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Introduction 


The general objective of this project has been to advance our theoretical understanding of 
Io's atmosphere and how various atomic and molecular species are lost from this atmosphere and 
are distributed in the circumplanetary environment of Jupiter. This grant has provided support 
for the activities of Dr. Michael Combi at the University of Michigan to serve as a small part in 
collaboration with a larger project awarded to Atmospheric & Environmental Research, Inc., 
with primary principal investigator Dr. William H. Smyth. Dr. Combi is the Principal 
Investigator and Project Manager for the Michigan grant NAG5-6187. This Michigan grant has 
provided for a continuation of a collaboration between Drs. Smyth and Combi in related efforts 
beginning in 1981, and with the object to develop and apply sophisticated theoretical models to 
interpret and to relate a number of new and exciting observations for the atmospheric gases of 
the satellite. The ability to interpret and then to relate through the theoretical fabric a number of 
these otherwise independent observations are a central strength of this program. This 
comprehensive approach provides a collective power, extracting more from the sum of the parts 
and seeing beyond various limitations that are inherent in any one observation. Although the 
approach is designed to unify, the program is divided into well-defined studies for the likely 
dominant atmospheric gases involving species of the S0 2 family (SO z , SO, 0 2 , S and O) and for 

the trace atmospheric gas atomic sodium and a likely escaping molecular ion NaX + (where NaX 
is the atmospheric molecule and X represents one or more atoms). 

Program of Work during the First Two Years 

Work during the first two years of this program involved three major efforts. The first 
was adapting our model (Smyth and Combi 1991) for the extended neutral sodium corona 
(zenocorona or magnetonebula) for use in calculating the distribution and escape of atomic 
oxygen and sulfur due to the fast escape from Io from charge exchange (-50-100 km/s) with fast 
torus ions and direct ejection from knock-on processes (-15-50 km/s). The second was using the 
near sodium corona model (Smyth and Combi 1988a&b) for understanding some of the new the 




eclipse and emission images of sodium by Schneider et al. (1990). We published a rather 
extensive paper on these results (Smyth and Combi 1997) during the first year, and a reprint of 
this paper is attached to this report. The third involved construction of our model for the 
molecular-ion source of sodium using our newly developed particle trajectory calculations for 
probing the physics of the nature of the molecular-ion source (Schneider et al. 1991). Our 
portion of this work was completed, and has already been reported in separate annual reports to 
NASA. 

Program of Work for the Third Year 

Our role in this project has been determined largely in response to the direction of Dr. 
Smyth at AER, Inc. who is the main Principal Investigator of the overall project. Work for this 
third and final year of the project involved a shift in our planned role in the program. The work, 
however, still falls within the specific area of studies for the loss of atomic and molecular species 
from Io. During the last two years Dr. Smyth and co-workers have developed a two-dimensional 
axisymmetric model for the detailed distribution of neutral species of the S0 2 family (namely 
S0 2 , SO, 0 2 , S and O) in Io’s atmosphere. The purpose of this is to provide a global picture of 
the distribution of atmospheric constituents from the surface, through the extended corona, and 
finally out to the escaping neutral gas in the circumjovian environment. More specifically, we 
want to put a “first-principles” physical basis to the results of our previous work (Smyth and 
Combi 1997) which found a particular form for the velocity distribution of atomic sodium was 
required to produce the observed radial (altitude) distribution of sodium in Io’s corora. 

In order to calculate the distribution of neutral species, it is necessary to characterize the 
full plasma environment including ion density, temperature and velocity. Plasma conditions are 
needed to characterize neutral species ionization rates and ion-neutral heating through elastic and 
charge exchange collisions. At his direction we supplied model results using our new three- 
dimensional MHD model for the interaction of Io with the corotating Jovian plasma torus, 
developed under a different grant (Combi et al. 1998). This included the massive set of 




numerical results behind those shown in our paper, as well as comparable results from new 
model runs which have not yet been published. As part of this work, in addition to simply 
providing the model results, we have been consulting with him in their appropriate interpretation 
and use in the context of his neutral atmosphere model, especially as they relate to the lower 
boundary conditions. Because his grant year lags ours by several months, it is not yet possible to 
show specific results of this work here, but we expect new results in the coming months. As an 
example of the type of MHD results we provided, we attach a reprint of our Io MHD paper 
(Combi et al. 1998) published last year. 
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A data set composed of different groundbased observations 
for lo’s sodium corona and spatially extended sodium cloud 
and covering the spatial range from lo’s nominal exobase of 

I. 4 satellite radii to east-west distances from Io of ± 100 satellite 
radii (R lo ) is used to investigate the velocity distribution of 
sodium at the exobase. The data set is composed of the novel 
1985 eclipse measurements of Schneider et ai (1991, Astrophys. 

J. 368, 298-315) acquired from —1.4 to —10 R| (P , the 1985 
east-west emission data of Schneider et al. acquired from —4 
to —40 R\ 0 , and sodium cloud image data acquired near lo’s 
orbital plane from —10 to —100 Ri 0 by a number of different 
observers in the 1976 to 1983 time frame. A one-dimensional 
east-west profile that contains Io is constructed from the data 
set and is analyzed using the sodium cloud model of Smyth 
and Combi (1988, Astrophys. J. Supp. 66, 397-411; 1988, 
Astrophys. J. 328, 888-918). When the directional feature in 
the trailing cloud is either north or south of this east-west line 
(i.e., not at the null condition), an isotropic modified [incom- 
plete ( a = 7/3) collisional cascade] sputtering flux speed distri- 
bution at the satellite exobase with a peak at 0.5 km sec -1 
provides an excellent fit to the data set for a sodium source of 
1.7 x 10 26 atoms sec -1 . In particular, the model calculation 
reproduces (1) the essentially symmetric column density distri- 
butions exhibited by the eclipse measurements about Io within 
the Lagrange sphere radius (5.85 R ]f} , i.e., the gravitational grasp 
of the satellite), (2) the change in the slope of the column 
density observed just beyond the Lagrange sphere radius in the 
east-west profile of the forward cloud, but not in the trailing 
cloud, and (3) the distinctly different east-west brightness pro- 
files exhibited by the forward and trailing clouds in the emission 
data at the more distant (— ±20-100 R ln ) portions of the cloud. 
In contrast, the speed dispersion at the exobase for either an 
isotropic Maxwell-Boltzmann flux speed distribution or an iso- 
tropic classical (a = 3) sputtering flux speed distribution (which 
has a higher velocity-tail population than the Maxwell- 
Boltzmann, but not as high as the incomplete collisional cascade 


sputtering distribution) is shown to be inadequate to fit the 
data set. To fit the enhanced trailing east-west profile observed 
when the directional feature is at the null condition, an addi- 
tional enhanced high-speed ( — 15-20 km sec -1 ) sodium popula- 
tion is required which is nonisotropically ejected from the satel- 
lite exobase so as to preferentially populate the trailing cloud. 
The need for such a nonisotropic high-speed population of 
sodium has also been recognized in the earlier modeling analysis 
of the directional features (Pilcher et al ., 1984, Astrophys . J. 
287, 427-444), in the more recent lower-velocity component 
required in modeling the sodium zenocorona (Smyth and 
Combi, 1991, J. Geophys. Res. 96, 22711-22727; Flynn et al ., 
1992, Icarus 99, 115-130), and in the very recent modeling of 
the directional feature reported by Wilson and Schneider ( 1995, 
Bull. Am. Astron. Soc. 27, 1154). A complete sodium source 
rate speed distribution function at lo’s exobase from 0-100 km 
sec -1 is then constructed by combining the isotropic modified 
[incomplete (a = 7/3) collisional cascade] sputtering flux speed 
distribution, the nonisotropic directional feature (lower-velocity 
zenocorona) source (—15-20 km sec 1 ) , and the higher-speed 
(—20-100 km sec -1 ) charge-exchange source required to simu- 
late the sodium zenocorona far from Jupiter, o mi Academic Press 


1. INTRODUCTION 

Atomic sodium in the Jupiter system originating from 
a satellite source at Io has been observed in the D : (5889.95 
A) and D { (5895.92 A) emission lines during the past 25 
yr from groundbased facilities. Using an observing slit, 
the sodium emission which is excited by solar resonance 
scattering was first discovered in 1972 by Brown (1974) 
very near lo, where its intensity is brightest [—many tens 
of kiloRaleighs (kR)] and where the sodium density is 
dominated by low-speed (—2 km sec 1 or less) ballistic 
atom orbits in the satellite “corona." By occultation of the 
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bright region near Io, image observations (Murcray 1978; 
Murcray and Goody 1978; Matson et al. 1978) were first 
acquired in 1976 and 1977 for fainter (— few to —0.5 kR) 
sodium more distant from Io but still near its circular orbit 
(radius of 5.9 Jupiter radii) about the planet and revealed 
the presence of a predominant “forward cloud 1 ’ and a less 
spatially extensive “trailing cloud” that moved with the 
satellite. This sodium has been characterized primarily by 
a source of low-speed (—2.6-4 km sec -1 ) atoms that have 
sufficient energy to just escape from Io with an excess 
velocity of only —1 km sec" 1 (or so) and thereby remain 
gravitationally bound to Jupiter fairly near the satellite 
orbit. Additional observations (Pilcher et al. 1984; Gold- 
berg et al. 1984) of even fainter (—1 to 0.2 kR) sodium in 
the early 1980’s revealed a “directional feature” attached 
to Io in the trailing cloud that oscillated north and south 
about the satellite plane with a phase and period deter- 
mined by the Io System III longitude angle. This sodium 
source was characterized by atoms with speeds —20 km 
sec -1 ejected nonisotropically from the satellite so as to 
populate the trailing cloud and the circumplanetary space 
at larger radial distances beyond Io’s orbit. From earlier 
slit measurements in 1974 (Trafton and Macy 1978), fainter 
(—30 R) sodium emissions well beyond Io’s orbit had been 
observed at a radial distance of —60 planetary radii, while 
from more recent images (Mendillo et al. 1990), very faint 
(—1 R) sodium emissions were observed at radial distances 
of —400-500 planetary radii. Sodium at these larger radial 
distances is called the “magneto-nebula” or “sodium zeno- 
corona” and is thought to be populated primarily by a 
nonisotropic charge-exchange source of high speed (—15- 
100 km sec l ) atoms at Io with velocity skewed in the 
forward direction of corotational plasma motion past the 
satellite, and secondarily by a narrow forward sodium jet 
produced by a spatially distributed molecular ion source 
(Wilson and Schneider 1994). Most of this sodium escapes 
the Jupiter system, forms a sodium pause in the sunward 
direction at —2300 planetary radii because of solar radia- 
tion acceleration, and is eventually lost to the solar wind 
by photoionization (Smyth and Combi 1991). 

The observations of sodium emissions on many different 
spatial scales in the Jupiter system thus indicate that its 
atomic source at Io’s exobase must have a wide dispersion 
of speeds. Modeling of these observations has in the past 
been mostly undertaken separately for only one of these 
spatial regions at a time. Although the higher velocity 
dispersions for the sodium zenocorona may be reasonably 
well understood because of its large spatial structure and 
the lack of any significant sodium lifetime impact of the 
magnetosphere, a consistent source for the slower sodium 
in lo's corona and in the forward and trailing clouds near 
its orbit has not been established. The recent determination 
of the sodium spatial profile in the Io corona obtained 
from the groundbased eclipse data of Schneider (1988; 


Schneider et al. 1987, 1991) coupled with earlier emission 
observations, however, now provides a viable observa- 
tional base from which it is possible to pursue the nature 
of this slower sodium. The investigation of a consistent 
exobasc sodium source for Io’s corona and the forward 
and trailing clouds near its orbit is therefore undertaken 
in this paper. A consistent flux speed distribution at the 
exobase is determined, and the corresponding sky-plane 
spatial distribution of sodium near Io is presented. Sodium 
source information obtained from previous modeling anal- 
ysis of Io’s corona and the forward and trailing clouds is 
first summarized in Section 2. The observational data base 
to be investigated in this paper is presented in Section 3. 
Modeling of an east-west spatial profile determined from 
this observational data base is undertaken in Section 4. 
Discussion and conclusions are presented in Section 5. 

2. EARLIER SODIUM MODELING 

The major modeling analysis studies for the spatial distri- 
bution of sodium near Io and its orbit are summarized in 
Table I. The summary is divided into three observed spatial 
regions: (1) the Io corona located within the satellite La- 
grange sphere (average radius of 5.81 R[ t) or —3 arcsec), 
(2) the neutral cloud located beyond the Lagrange sphere 
and near Io’s orbit, and (3) the north-south oscillating 
directional feature, observed to trail Io in its orbit. Model- 
ing analysis studies for the Io corona are further subdivided 
into early observations of the average intensity in an 8 X 
3 arcsec slit centered on Io reported by Bergstralh et al. 
(1975, 1977) that indicated an east-west intensity asymme- 
try of —1.25 and later observations for one-dimensional 
column density profiles within the Lagrange sphere re- 
ported by Schneider et al. (1987, 1991). 

2.1. Corona: East-West Intensity Asymmetry 

In Table I, the early studies of Smyth (1983) for sodium 
atoms ejected monoenergetically from Io’s exobase estab- 
lished that small scale structures in the D-line intensity 
profile observed as a function of the Io geocentric phase 
angle (Bergstralh etal. 1975, 1977) could arise from modu- 
lation of the atoms’ escape rate from Io caused by the 
action of solar radiation acceleration in the D-lines. These 
modulations occur primarily for exobase speeds near 2.0 
and 2.1 km sec \ which arc near the escape-speed thresh- 
old of the Lagrange sphere. Later studies of Smyth and 
Combi (1987a) showed that the main reason for the east- 
west intensity asymmetry was, however, an east-west elec- 
tric field which altered the plasma properties at Io’s orbit 
so as to increase the sodium lifetime and hence sodium 
abundance when Io was preferentially east of Jupiter. More 
complex modeling studies of Smyth and Combi (1988b) 
constrained the flux velocity dispersion for sodium at Io’s 
exobase by simultaneously fitting the average east-west 



TABLE I 

Summary of Modeling Studies for the Spatial Distribution of Sodium Near Io and Its Orbit 
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intensity asymmetry and also the general spatial morphol- 
ogy of the forward sodium cloud, located on a much larger 
spatial scale well beyond the Lagrange sphere. These stud- 
ies showed that the sodium ejection speed at the exobasc 
required to fit the east-west intensity asymmetry is double- 
valued, having a lower value of <1 km sec 1 and a higher 
value in the range 2.6-3.65 km sec -1 . For a Maxwell- 
Boltzmann flux distribution, the lower and higher most 
probable speed values were 0.71 km sec -1 ( T = 460 K) 
and 3.65 km sec -1 ( T = 12,300 K). Neither distribution 
was, however, suitable for properly populating the forward 
cloud. The lower value produces essentially only ballistic 
atom orbits which could not populate the forward cloud, 
while the higher value was significantly larger than the 
nominal —2.6 km sec -1 characteristic monoenergetic veloc- 
ity required to reproduce the proper spatial morphology 
of the forward cloud as a function of the Io geocentric 
phase angle. For a Maxwell-Boltzmann flux distribution 
with a more nominal thermal exobase temperature in the 
range —1000-2000 K, the calculated east-west intensity 
ratio was much higher than the observed value, with the 
atoms still contributing primarily to the corona density 
and again far too deficient in energy to contribute any 
significant sodium to the forward cloud. For the preferred 
(a; = 7/3) modified-sputtering distribution of Smyth and 
Combi (1988b) with source strength —2 X 10 26 atoms sec -1 , 
the lower and higher speed values were <0.5 and —2.9 km 
sec -1 , respectively, with the latter value being preferred 
because of its closer proximity to the —2.6 km sec - 1 charac- 
teristic velocity for the forward cloud. Interestingly, how- 
ever, it is actually the lower value that will be shown in 
this paper to reproduce the correct spatial profile for so- 
dium both within the Lagrange sphere and beyond in the 
more distant neutral cloud. 

2.2. Corona: Column Density Profile 

In Table I, modeling studies of Smyth and Combi 
( 1987b, c) determined that typical forward cloud brightness 
data for the sodium cloud could be properly simulated 
well beyond the Lagrange sphere radius of —5.81 R Ul by 
a sodium source of —1 X 10 26 atoms sec -1 ejected monoen- 
ergctically from Io’s exobase with a characteristic velocity 
of —2.6 km sec -1 . They also established that this same 
sodium source reproduced the column density profile of 
Schneider et al. (1987) within the Lagrange sphere down 
to a radius of —3.5 R lo . For a radius smaller than —3.5 /? lo , 
the calculated profile was lower than the observed profile, 
indicating that lower (ballistic) velocity components are 
required, in addition, as part of a more realistic flux velocity 
dispersion. A similar behavior for the simulated column 
density profile, with an even more dramatic departure from 
the observed profile both inside and outside the Lagrange 
sphere, was also later shown by a model calculation of Ip 


(1990), who assumed an exobase speed of 3 km sec 5 but 
did not include the gravity of Jupiter so as to properly 
include the near zero escape speed conditions for sodium 
at the Lagrange sphere. Adopting for sodium atoms at the 
exobase a simple (i.e., binding-energy velocity v b = 0) 
classical sputtering energy distribution with a low energy 
cut-off and also excluding Jupiter’s gravity, McGrath 
(1988) modeled the column density within the Lagrange 
sphere and produced a profile with a slope slightly less 
steep than the observation for an infinite sodium lifetime 
and a slope somewhat steeper than the observation for a 
sodium lifetime of 3 hr. Alternatively adopting a Maxwell- 
Boltzmann flux distribution, assuming an infinite sodium 
lifetime, and similarly excluding Jupiter’s gravity. Summers 
etai (1989) and Schneider etai (1991) modeled the column 
density within the Lagrange sphere region and produced a 
profile that reasonably well matched the observed profile 
for an exobase temperature, respectively, of 1000 K based 
on the partial eclipse data set (Schneider et al. 1987) and of 
1500 K based on the complete eclipse data set (Schneider et 
ai 1991). Although these different flux velocity distributions 
reasonably fit the observations within the Lagrange sphere, 
it is clear from the earlier studies of Smyth and Combi 
(1988b) that these Maxwell-Boltzmann distributions are 
energetically deficient and inappropriate for populating the 
neutral cloud and, furthermore, that the more energetically 
promising sputteringdistribution cannot be investigated ad- 
equately near or beyond the Lagrange sphere radius without 
properly including the gravity of Jupiter, solar radiation ac- 
celeration, and the spacetime variable sodium lifetime in the 
plasma torus. This study will be undertaken in Section 4. 

2.3. Sodium Cloud 

The early studies in Table I for the sodium cloud were 
general in nature, probing its poorly documented spatial 
and angular extent about the planet. Based upon the solar 
resonance scattering excitation mechanism for sodium 
(Bergstralh et al. 1975) and limited angular extent data 
determined by slit-averaged intensity data, Carlson et al. 
(1975) undertook monoenergetic (3.5 km sec ! ) model cal- 
culations and estimated that the sodium cloud lifetime 
(assumed to be spatially uniform) was likely determined 
by electron impact ionization by the (then very poorly 
characterized) plasma in the planetary magnetosphere. 
This general picture for the cloud was confirmed by more 
extensive model calculations performed by Fang et al. 
(1976) and Smyth and McFJroy (1977), the latter of which 
explored the time evolution and two-dimensional nature 
of the cloud for exobase velocities near the Io Lagrange 
escape speed. The acquisition of sodium cloud images in 
late 1976 and early 1977 brought this subject into dramatic 
focus. For a classical sputtering flux distribution that 
peaked at 4 km sec' 1 , Matson et al. (1978) successfully 
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modeled a one-dimensional east-west brightness profile 
(derived from a cloud image) which extended from lo in 
the forward cloud to -80 /? lo and in the trailing cloud to 
—40 Rio but which excluded sodium emission within Io’s 
corona. The analysis (Smyth and McElroy 1978) of a much 
larger sodium cloud image data set (Murcray 1978) also 
indicated that the forward cloud could be characterized 
by an exobase ejection speed of —2.6 km sec 1 and that 
its changing intensity pattern could be understood as the 
changing viewing perspective of an approximately steady 
state cloud on the sky plane as lo moved on its orbit around 
Jupiter. The observed predominance of the forward so- 
dium cloud over the trailing cloud was accomplished in all 
these models by limiting the exobase source area to a 
hemisphere (see Table I) and by limiting the assumed 
spatially uniform lifetime so as to dynamically select so- 
dium atom orbits that would primarily populate the for- 
ward cloud. Additional modeling studies by Macy and 
Trafton (1980) of the radial and vertical cloud structure 
on a larger spatial scale indicated source dispersion speeds 
at least up to 13 km sec 1 were, however, required to 
explain a variety of other observations. Additional model 
studies (Smyth 1979, 1983) showed that the newly discov- 
ered east-west orbital asymmetry of the sodium cloud 
(Goldberg et al. 1978) was not source related but was due 
to the perturbing action of solar radiation acceleration 
on the sodium atom orbits. Adopting a one-dimensional 
radially dependent sodium lifetime in the plasma torus 
based upon limited Voyager spacecraft data and an asym- 
metric exobase source for a classical sputtering distribution 
with a peak velocity at 4 km sec 1 , Goldberg et al. (1980) 
successfully modeled a one-dimensional east-west bright- 
ness profile acquired during the Voyager 1 encounter for 
distances extending from lo in the forward cloud to —80 
R lo and in the trailing cloud to —30 R} v . Later modeling 
by Smyth and Combi (1988b) using a more accurate two- 
dimensional space- and time-dependent sodium lifetime 
in the plasma torus and an isotropic (or near isotropic) 
exobase sodium source of —2 X 10 26 atoms sec 1 demon- 
strated that the predominant forward cloud was caused by 
the highly radially dependent sink for sodium in the plasma 
torus and not by a nonisotropic source. The deduced char- 
acteristic or most probable exobase speed for the more 
definitive modeling of the forward sodium cloud above is, 
therefore, in the range —2.6 to 4 km sec 1 and is much 
larger than required to characterize the sodium column 
density profile in Ios corona. A new flux speed distribution 
is therefore needed for consistency and is determined in 
Section 4. 

2. 4. Directional Feature 

In Table I, observations acquired in 1980 and 1981 by 
Pilcher et al. (1984) for weaker D-line emissions in the 


trailing portion of the sodium cloud allowed them to dis- 
cover an elongated feature in the brightness distribution 
that on the sky plane was directed away from Jupiter and 
was inclined sometimes to the north and sometimes to 
the south of the satellite’s orbital plane. The north-south 
direction of the feature was shown to be correlated with 
Io’s magnetic longitude and suggested a formation mecha- 
nism involving the oscillating plasma torus. Modeling anal- 
ysis by Pilcher et al. indicated that the feature resulted 
from a high-velocity (—20 km sec' 1 ) sodium source that 
was at near right angles to Io’s orbital motion with a source 
strength required on the outer satellite hemisphere of 
— 1 x 10 26 atoms sec" 1 . This peculiar directionality of the 
source was investigated by Sieveka and Johnson (1984), 
who concluded that it was likely produced by direct colli- 
sional ejection of neutral sodium from the exosphere by 
the corotating plasma flow past lo. Modeling of the sodium 
zenocorona (Smyth and Combi 1991; Flynn et al. 1992) 
showed that it was consistent with a two-component exo- 
base source: a similar high-velocity (—20 km sec' 1 ) sodium 
source of —1 X 10 26 atoms sec" 1 for the spatial distribution 
nearer the planet and an even higher-velocity (—57 km 
sec' 1 ) sodium source of —2 x 10 26 atoms sec 1 for the 
spatial distribution further from the planet. Both source 
components, however, were based on ion-neutral charge 
exchange processes in Io’s exosphere and were hence com- 
posed of a speed tangential to Io’s orbit at Io’s position 
plus an isotropic Maxwell-Boltzmann distribution with a 
most probable speed of about one-third of the tangential 
speed in the Jupiter frame (i.e., one-third of —37 and —74 
km sec' 1 , respectively). In Table I, the lower component 
is therefore symbolically denoted in the lo frame by 20± 12 
km sec' 1 . Recent modeling of the directional features has 
also been reported by Wilson and Schneider (1995), who 
used a similar lower component source denoted in Table 
I by 20 ± 10-20 km sec' 1 , where the isotropic portion of 
their source may be variable in magnitude. A spatially 
distributed molecular-ion source has also been used to 
model the sodium zenocorona (Flynn 1993), but appears 
to be of secondary importance due to its smaller source 
rate and erratic presence. 

3. OBSERVATIONAL DATA BASE FOR MODELING 

To describe the entire spatial distribution of sodium in 
Io’s corona and beyond in the extended neutral clouds, 
three different types of sodium observations obtained on 
very different spatial scales are combined. For Iocentric 
distances, the combined data set is composed of the novel 
1985 eclipse measurements of Schneider et al. (1991) ac- 
quired from —1.4 to —10 R\ 0 , the 1985 east-west emission 
data of Schneider et al. (1991) acquired from -4 to -40 
R lo , and sodium cloud image data acquired near Io’s orbital 
plane from -10 to -100 R io by a number of different 


a 
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TABLE II 


lo Eclipse and Emission Observations for 1985 


Date 

UT 

■Midpoint 

lo Geocentric 
Phase Angle 
Range 
(deg) 

lo System HI 
Longitude 
Range 
(deg) 

Spectrum ID 

Type of 
Observalion 

Eclipse Emission 

Dominant 

Spatial 

Profile 

Directional 

Feature 

Orientation 

Enhance 

Doppler 

Signature 

D; Emission Profile Power Law Fit * 
Exponent Amplitude (kR) 

East West East West 

Profile Profile Profile Profile 

August 27 

0714 

61.4 ± 

0.7 

29 9 t 2.3 

85gl88 

l 

forward 

null 

trai ling/forward 

1 67 (F) 

1.57 (T) 

191 

124 


072(1 

62.3 ± 11.0 

32 9 t 35.8 

a 

1 









0830 

72.2 ± 

1.1 

65.1 ± 3.4 

85gl96 

2 

symmetric 

south 

no 

1 85 (S) 

1 80 <S) 

169 

142 

September 1 3 

0641 

276 6 ± 

0.7 

194 7+ 2 3 

85h032 

3 

trailing 

null/north 

trailing 





September 14 

0245 

87.7 ± 

07 

31.4 ± 2.3 

8 5 h 1 0 2 

4 

trailing 

null 

trailing 

1 231S/T) 

1 80 (S/F) 

89 

188 


0326 

93.5 t 

4.2 

50.5+13.7 

b 

2 









0416 

100 6 + 

0.4 

73.6+ 14 

85hl 13 

5 

8 

south 

no 

1 27 (T) 


66 


September 15 

0316 

294.7 ± 

0.7 

353.7+ 3.0 

85h 1 52 

6 

trailing 

null 

(railing 






0500 

309.4 ± 

11.5 

41.7 + 37.4 

c 

3 


— 






September 21 

0604 

100.5 + 

1.6 

112.5+ 5.2 

d 

4 








September 23 

0230 

1 17.2 ± 

0.7 

267.61 2 3 

85h433 

7 

forward 

north 

no 

2 05 (T) 

1 57(F) 

342 

135 


0301 

121.6 + 

0.4 

281.9+ 12 

85h436 

8 

forward 

north 

no 

1.96 (T) 

1 54(F) 

283 

138 


0356 

129.3 ± 

20 

307.3+ 6.7 

e 

5 









0534 

143.1 ± 

07 

352.9 t 2.3 

85h457 

9 

forward 

null 

no 

2 .16 (T) 

1 64 (F) 

374 

165 


a. Eclipse 1 B5g179, 85g181, 85g185, 85g 1 Sa 05g192. 85g193. 05g196 

b Eclipse 2 B5h103. 85h1 04. 856105, 85h106, 85MC7. 85MO0, 056109. 856110, 856112 

C Eclipse 3 85h1 53, 856154 856155, 856157, 856159. 8Sh162. 85h163 
d Eclipse 4 856207, 856280 856289, 85h290. 856291 656292. 856293. 85h294. 856295 

e. Eclipse 5 856441. 856442. 85h443, 056444 056445. B56446. 056447. 856448 85h449. 856450 

f. Profile points inside of 4 R (o arc excluded, power law fu A r'P, where A is the amplitude, p is the exponent, and r is in units of Rj 0 ; F - forward cloud; S = symmetric turning point. T - trailing cloud 

g. Nol sufficient data west of lo lo compare spatial profits (sec Table III) . 

(Note one lower-bound data point from eclipse 4 at a distance from the center oflo of 1.17 Rj 0 is excluded in ihe analysis since it is well within the nominal exobase radius of 14 R| 0 ) 


observers in the 1976 to 1983 time frame. For modeling 
purposes in Section 4, a one-dimensional east-west profile 
centered on lo is constructed from this data set. From the 
data of Schneider et ai (1991), five higher quality eclipse 
profiles and nine higher quality emission profiles have been 
selected, and their observational dates, times, lo angular 
parameters, spectral ID numbers, and the numbering of 
these interleaved observations as adopted in this paper 
are summarized in Table II. The D 2 brightnesses for the 
emission profiles in Table II, previously published only in 
a graphical format, are given numerically in Table III as 
provided by Schneider (1990, 1995, both private communi- 
cations). For the sodium cloud data, fourteen images in the 
D 2 emission line acquired in the 1976-1983 time interval 
(Murcray 1978; Murcray and Goody 1978; Matson et ai 
1978; Goldberg et ai 1980, 1984; Morgan 1984, private 
communication) have been selected, with values for the 
east-west D 2 brightness profiles of the forward and trailing 
clouds extracted and summarized in Table IV. 

The eclipse observations provide the most accurate in- 
formation in the radial interval from —1.4 to 6 /? Io for the 
atomic sodium column density profile in the corona within 
the Lagrange sphere of lo ( i.e., a radius of 5.81 R io ) and 
yield an essentially symmetric column density profile about 
lo with a power law lit N (1.4 < r < 5.85) = 2.55 X 
If)! 2 /* 248 , where N is in units of atoms cm 2 and r is the 
distance from the Io’s center in units of R lo . This power 
law fit, however, undercuts the eclipse data beyond the 
Lagrange radius, and this reduced slope will later be seen 
to be caused by the dominant planetary gravitational field 


beyond Io’s Lagrange radius. The emission observations 
provide accurate information for the sodium D 2 brightness 
from just within the Lagrange sphere outward into the 
nearer portion of the sodium cloud (i.e., —4-40 R l0 from 
Io’s center) with a power law fit / D ^(r^4) = 101 r“ 145 , 
where I D ^ is in units of kiloRayleighs (kR) and where the 
D 2 brightness of —100 kR as r approaches Io’s surface is 
consistent (see Brown and Yung 1976) with the maximum 
sodium column density of —1 X 10 12 atoms cm" 2 deduced 
from the eclipse data. The brightnesses for the different 
observed profiles in Table III, however, vary by a factor 
of —3 to almost 5 at the same distance from lo but have 
error bars that are no larger than ±30%, suggesting that 
the large variation is real and likely correlated with lo 
geocentric phase angle, lo System III longitude, and the 
east-west asymmetry in the plasma torus, as is the case 
for sodium cloud image data. The closer spatial regions 
covered by the eclipse data and the near lo emission data 
are, however, masked in the cloud images by a circular 
(or nearly circular) occulting mask of —10 R Ui in radius 
centered on lo. The sodium cloud image data in Table IV, 
as illustrated in Fig. 1, hence contribute to the east-west 
brightness profiles for Iocentric distances from —10 to 
-100 R la . 

The structure of the sodium cloud emission brightness 
on the sky plane has been historically divided into a for- 
ward cloud, so-called because it appears ahead of the satel- 
lite in its orbit (i.e., in Fig. I located right (west) of lo in 
image A and left (east) of lo in images B and C), and 
a corresponding trailing cloud that appears behind the 
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TABLE III 

Emission Data for 1985 


Radial 
Distance 
From Io 
(satellite radii) 

D2 Intensityt(kR) 

Emission 1 
61.4° 

27 August 
(85gl 88) 

Emission 2 
72.2° 

27 August 
(85gl96) 

Emission 4 
87.7° 

14 September 
(85h 102) 

Emission 5 
100.6° 

14 September 
(85h 1 13) 

Emission 7 
117.2° 

23 September 
(85h433) 

Emission 8 
121.6° 

23 September 
(85h436) 

Emission 9 
143.1° 

23 September 
(85H457) 

Emission 3* 
276.6° 

13 September 
(85h032) 

Emission 6 
294.7° 

15 September 
(85hl52) 


-29.92 





0.69 ±0.17 

0.93 10.18 

0.6610.16 




-29.65 









1.1210.19 


-24.48 








1.4010.21 



-21.76 





1.1910.19 

0.8410.17 

0.95+0.18 




-16.32 





1.3810.20 

2.1010.26 


2.22 ±0.27 



-15.78 







1.7410.23 


2.1510.26 


-14.96 




1.91 ±0.24 







-10.88 

3.16 ± 0.35 




2.9010.33 

3.2910.36 

3.3810.37 

4.40 ± 0.46 


T 

-9.52 



3.1710.35 







West 

-6.80 

5.21 ±0.54 

4.50 ±0.47 

6.41 ±0.66 


7.4910.76 

7.76 1 0.79 


11.541 1.16 



-4.08 

14,65 ± 1.47 

1 1.29 ± 1.14 

* 14.71 ± 1.48 




7.54 ± 0.77 

28.6912.87 



-3.81 









18.75 12.63 


-2.18 

24.31 + 3.40 

18.39 ±2.57 

24,1213.38 


25.61 13.59 

25.01 ±3.50 

21.0012.94 

45.4716.37 



-1.09 

27.46 ± 3.84 

19.99 + 2.80 

27.1413.80 


29.95 14.19 

28.72 + 4.02 

23.89 ±3.34 

48.61 ± 6.81 



-0.82 









24.85 1 3 48 

Io 

0.00 

27.44 ± 3.84 

19.19 + 2.69 

27. 14 ± 3.80 


27.9013.91 

30.1614.22 

23.97 ± 3.36 

45.5616.38 



0.27 









25 27 1 3.54 


1.09 

25.87 ± 3.62 

19.35 ±2.71 

26.65 ± 3.73 

23.03 1 3.22 

26.31 ±3.68 

26.39 1 3.69 

22.84 ± 3.20 

39.45 + 5.52 



1.36 









22.901 3.21 


2.18 

25.19 ± 3.53 

16.87 ±2.36 

25.68 ± 3.60 

19.60 + 2.74 

22.42 ±3.14 

23.91 ±3.35 

20.21 ±2.83 

33.6014.70 



3.26 









11.53 1 1.61 


4.08 

17.60 ± t .77 

1 1.92 ± 1.20 

18. 10 ± 1.82 

11.871 1.20 

15.541 1.56 

14.831 1.49 

12.52 ± 1.26 

20.0612.01 



5.98 









4.7910.50 

East 

6.80 

8.31 ±0.84 

5.52 ±0.57 

8.1710.83 

5.5210.57 

7.0910.72 

7.1010.73 

6.0210.62 

7.5210.77 


i 

10.06 









1.6610.22 


10.88 

3.40 + 0.37 

1.94 ±0.25 

3.8410.41 

3.21 ±0.35 

2.47 10.29 

2.47 1 0.29 

2.2610.27 

2.73 ±0.31 



15.50 









0.77 ±0.17 


16.32 



2.8410.32 

1.81 ±0.24 

1.03 ±0.18 

1.11 ±0.19 

0.8210.17 

1.3910.20 



21.76 



2.25 10.27 


0.68 10.16 

0.7510.17 

0.5210.16 




23.66 









0.5310.16 


24,48 




1.2410.19 




1.0610.18 



27.20 



1.3910.21 








35.36 



1.2710.20 





1.0010.18 



46.24 








0.83 ±0.17 



^Nonuniform spatial coverage occurs because of different distance intervals adopted to obtain good average brightness values (given different signal to noise ratios) and because 
of signal drop-oul associated with constraints imposed on positioning the slit profile on the CCD detector during interleaved eclipse and emission measurements. 

* Calibration uncertain 


TABLE IV 


East-West D 2 Brightness Profiles for Sodium Cloud Image Data 


UT Date ITT Time 

Image ID 
Number 

Io Geocentric 
Phase Angle 

(deg) 

lo System III 
Longitude 
(deg) 


East-West Distance from lo for Specified D 2 Bnghtness Level (satellite radii) 


Image Data Set Reference: Murcray ( 1 978) 



0.5 kR 


1 OkR 


1 5 kR 


2 OkR 






Forward 

Trailing 

Forward 

Trailing 

Forward 

Trailing 

Forward 

Trailing 

1976 Nov 16 0806 

ES 328B 

256 

262 



51 

29 

46 

24 

28 

14 

l 977 Jan 27 0024 

ES 369A 

86 

193 

60 

38 

50 

22 

44 

19 

32 

16 

0217 

ES 370D 

102 

245 

>83 

30 

56 

25 

46 

21 

24 

19 

Image Data Set Reference: Goldbere( 1 988 1 ) 



0.2 kR 


0.5 kR 


TO kR 


2.0 kR 






Forward 

Trailing 

Forward 

Trailing 

Forward 

Trailing 

Forward 

Trailing 

1981 May 5 0819 

SLR 418/31-33 

102 

300 

78 

69 

63 

40-66 

45 

37 

24 

21 

May 1 2 0848 

SIP 420/30 32 

91 

302 

66 

86 

41 

37 

35 

29 

20 

21 

May 1 3 0346 

SIP 42 1/2 1-23 

253 

108 

74-103 

73 

70 

41 

41 

27 

29 

20 

0555 

SIP 421/32-33 

271 

168 

124 

65 

112 

44 

51 

30 

26 

19 

June 6 0436 

SEP 424/10-12 

103 

300 

81 

75 

68 

71 

36 

26 

20 

23 

Image Data Set Reference. Morgan (1984+) 



0 3 kR 


0.6 kR 


0,9 kR 


1.8 kR 






Forward 

Trailing 

Forward 

Trailing 

Forward 

Trailing 

Forward 

Trailing 

1983 June 13 0714 

t 8492 

274 

230 


39 

58 

26 

37 

19 

20 


0722 

» 8494 

275 

233 


35 

52 

22 

40 

20 

20 


0729 

i 8496 

276 

237 

>93 

36 

52 

22 

37 

17 

22 


0827 

i 8501 

284 

264 

>93 

39 

61 

25 

42 

19 

23 


0949 

i 8509 

296 

301 

>93 

47 

63 

26 

44 

17 

24 


1010 

i 8511 

299 

311 

>93 

41 

63 

23 

44 

16 

19 



Umage observational data obtained by private communication 
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10 SODIUM CLOUD 1981 



FIG. 1. Io sodium cloud images. Three calibrated D 2 emission images of the lo sodium cloud from the JPL Table Mountain Data Set are shown 
to proper scale with Jupiter and Io’s orbit as viewed from earth in 1981 (Smyth and Goldberg 1993). The Io System III longitude and corresponding 
orientation of the trailing directional feature in image A are 247° and north, in image B are 104° and south, and in image C are 178° and only very 
slightly north. An east-west spatial scale of ±100 satellite radii about Io is shown for reference, with tick marks also located at ±50 satellite radii. 
Contour levels for the D 2 brightness, from outside to inside, are 0.2, 0.5, 1, 2, 4, 6, 8, and 10 kR. An occulting mask of —10 R\ a in radius is centered 
on Io so that brightness values within this distance are not accurate. 


satellite. The forward cloud changes in length, brightness, 
and east-west orientation relative to the satellite’s location 
as Io moves about Jupiter (i.e., as a function of the Io 
geocentric phase angle). This change in orientation is well 
documented (Murcray 1978; Murcray and Goody 1978; 
Goldberg etal 1984) and is due primarily to the projection 
upon the two-dimensional sky plane of a three-dimensional 
cloud (Smyth and McElroy 1978), slightly altered from 
mirror symmetry by solar radiation pressure (Smyth 1979, 
1983), which passes through an east symmetric turning 
point between about 65° and 85° to 90° and a west symmet- 
ric turning point at about 235° (Goldberg et al. 1984). A 
detailed examination of the emission data of Schneider et 
al. (1991), summarized in the last seven columns of Table 
II, shows that the forward and trailing profiles for the 
emission data are quite consistent with this known behavior 
of the Io sodium cloud images. In the trailing cloud, the 
time-dependent change in the north-south inclination of 
the fainter directional feature (see Fig. 1) has been shown 
to be correlated with the System III longitude of Io (Pilcher 
etal. 1984; Goldberg etal. 1984) with the directional feature 
changing from a south to north inclination (a first null 
point) at an Io System III longitude near 165° and changing 
from a north to south inclination (a second null point) for 


a rather poorly defined Io System III longitude somewhere 
between about 320° and 25°. When the directional feature 
is near the null location, as illustrated in Fig. 1, an increase 
in both the spatial extension and brightening of the trailing 
cloud along the east-west oriented (dashed) line is readily 
apparent. In addition, since the trailing cloud is associated 
with a high-speed Io sodium source (—15-20 km sec" 1 ), 
an increase in the Doppler width of the spectral line in 
the trailing cloud brightness along an east-west slit is also 
expected near the null location and is indeed observed in 
the emission data of Schneider et al. (1991), as indicated 
in Table II. 

Information for the sodium D 2 brightness profiles in 
both the forward and trailing clouds is presented in Fig. 
2. The forward and trailing cloud orientation depicted in 
Fig. 2 is chosen for Io near eastern elongation in order 
to facilitate the comparison with the 1985 emission data 
profiles, mostly acquired for Io east of Jupiter. The five 
emission profiles of Schneider et al. (1991) for Io east of 
Jupiter are shown by different symbol together with their 
power-law fits (/ Di (r — 4) - Ar~&) given in Table II. For 
the Io sodium cloud images, the extracted east-west D 2 
brightness profiles in Table IV for the forward and trailing 
clouds are shown by shaded areas, which represent appro- 
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TRAILING CLOUD 


FORWARD CLOUD 



FIG. 2. East and west brightness profiles for selected 1985 emission data and image cloud data. The spatial profiles both east and west of lo 
for the sodium D 2 emission brightness in units of kilorayleighs are shown as a function of the distance along the observing slit from the center of 
lo. Five emission observations identified by their satellite geocentric phase angles arc shown by the different symbols. These five profiles occur 
when lo is east of Jupiter and past the satellite phase angle where the forward cloud has its symmetric turning point so that the trailing cloud 
profiles are all to the east (left) of lo and the forward cloud profiles are all to the west (right) of lo. A power law fit to each profile is also shown. 
At larger distances from lo, an envelope for the east-west D 2 emission profile acquired from sodium image data is shown by the shaded area. For 
the trailing profile, the shaded area is divided into two parts: the lower area corresponding to sodium cloud data when the directional feature is 
oriented either north or south and the upper area corresponding to the directional feature oriented along the east-west direction (i.e.. the 
null condition). 


priatc bounds for the brightness profiles when lo is some- 
what near the elongation points of its orbit. For the trailing 
cloud, two different shaded areas are shown in Fig. 2 for 
the two different basic orientations of the directional fea- 
ture: (1) lower area, when the directional feature is inclined 
either north or south, and (2) upper area, when the direc- 
tional feature is at or near the null locations. As expected, 
the shaded area for the directional feature near the null 
locations is both brighter and less steep than the shaded 
area for the directional feature with either a significant 
north inclination or a significant south inclination. At larger 
distances from lo (>30 /? IO ), however, note that both 
shaded areas for the trailing sodium cloud are dimmer and 
more closely confined to lo than the shaded area for the 
forward cloud. In the forward cloud, all the emission 
brightness profiles are fairly tightly confined and have a 
slightly steeper slope than the shaded area due to increas- 
ing lo geometric plase angle. In the trailing cloud, the 
brightest and least steep of these profiles is for the emission 
4 (lo phase angle 87.7°) acquired for the directional feature 
at the null condition, while the next brightest profile is 
for the emission 5 (lo phase angle 100.6°), acquired ~~1.5 
hr later. 

4. ANALYSIS OF THE OBSERVATIONS 

Modeling analysis of the one-dimensional sodium distri- 
bution described in the previous section will now be under- 


taken. Collectively, the eclipse measurements for the co- 
rona near lo, the emission measurements that extend into 
the near sodium cloud, and the sodium cloud image derived 
profiles that reach to distances of ±100 7? Io , provide a set 
of spatially overlapping observations that will be used to 
study and constrain the initial velocity dispersion of the 
sodium source atoms at the exobase. In the modeling analy- 
sis, one-dimensional profiles are calculated using the nu- 
merical sodium cloud model of Smyth and Combi 
( 1988a, b), where the electron impact ionization sink for 
sodium is determined for a 7° tilted corotating plasma 
torus with an offset-dipole planetary magnetic field in the 
presence of a nominal (i.e., ~~2.8 mV m" 1 in Io’s frame) 
east-west electric field. A System III longitudinal asymme- 
try, although present in the torus ion emission, is not in- 
cluded but deferred to a later time when the electron de- 
pendence is available. 

To investigate the nature of the initial velocity dispersion 
of the sodium source, two different source flux speed distri- 
butions discussed earlier by Smyth and Combi (1988b; see 
their Appendix D) are considered: (1) a Maxwell- 
Boltzmann flux distribution and (2) a modified-sputtering 
flux distribution. The Maxwell-Boltzmann flux distribu- 
tion 4>( v\ T) is based on the Maxwell-Boltzmann velocity 
distribution and is defined as 

(I) 
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Speed (km/s) 

FIG. 3. Flux speed distribution functions for sodium at Io's exobase. Maxwell-Boltzmann flux speed distributions for sodium are shown for 
most probable speeds, v m , of 1.3 km sec -1 and 2.0 km sec -1 . Modified sputtering flux speed distributions are also shown for a = 3 and a most 
probable speed of 1.0 km sec 1 , for a 7/3 and a most probable speed of 0.5 km sec' 1 , and for a = 2 and a most probable speed of 0.4 km sec" 1 . 
All of the flux speed distributions are normalized to unit area under the curve, 


where v x = V2 kT/m is the most probable speed of the 
velocity distribution for an atom of mass ra. The Maxwell- 
Boltzmann flux distribution is proportional to the local 
velocity integrated flux referenced here to the satellite 
radius /? s not the exobase radius R E and depends upon 
one parameter, the exobase temperature T (or alterna- 
tively whic h determines both the most probable speed 
v m = V3kT/m and the speed dispersion of the flux distribu- 
tion. The modified-sputtering flux distribution 4>{v\ a, i> b , 
u M ) is proportional to the local velocity integrated flux 0 0 
and depends upon three parameters, an exponent a and 
two velocity parameters v h and p M , 

<f>(v;a,v b ,v M ) = 0o —— 1 — — (— ) 

\R e J v h D(a,v M fv b )\v b / 

V 1±^\ U2 ] 

w 2 m / J ’ 

where D(a , v M Iv h ) is a normalization constant (see Smyth 
and Combi 1988b). The exponent a primarily determines 
the dispersion of the distribution, which has a greater high- 
speed population as a decreases. The exponent a has a 
value of 3 for a classical sputtering distribution (i.e., a 


complete collisional cascade process) and a value of 7/3 
for a Thomas-Fermi modified-sputtering flux distribution 
(i.e., the limit of a single elastic collisional ejection process), 
where the latter distribution is based upon a Thomas- 
Fermi differential scattering cross section. The velocity 
parameter v b is related nonlinearly to the most probable 
speed v m of the flux speed distribution and primarily deter- 
mines v m (see Smyth and Combi, 1988b, Appendix D). The 
velocity parameter primarily determines the maximum 
speed for the flux distribution and depends upon the maxi- 
mum relative speed (and masses) of the plasma torus ion 
and sodium atom. For different values of their parameters, 
two Maxwell-Boltzmann flux distributions and three mod- 
ified-sputtering flux distributions are shown in Fig. 3 and 
will be utilized in the subsequent modeling analysis. 

In calculating the column density and the D 2 emission 
brightness in the numerical sodium cloud model, a smaller 
two-dimensional sky-plane grid centered on Io (±15 R lo ) 
is used to cover a spatial scale near the satellite more 
appropriate to the eclipse data while a much larger two- 
dimensional sky-plane grid centered on Io is used to cover 
a larger spatial scale more appropriate for the emission 
data and the sodium cloud image data. A one-dimensional 
profile for the eclipse data is obtained from the smaller 
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FIG. 4 . Model calculations for the lo eclipse data using a Maxwell- Boltzmann flux speed distribution. The atomic sodium column density 
profile near Io determined from the 1985 eclipse data by Schneider et al. (1991) is shown by the open circles. The model calculated column density 
profiles are shown by solid dots for the (cylindrically averaged) corona, by solid triangles for the forward cloud along the east-west slit direction, 
and by solid squares for the trailing cloud along the east-west slit direction. These column density profiles were calculated using the lo sodium 
cloud model of Smyth and Combi (1988b) for their case C description of the plasma torus and for an Io geocentric phase angle of 92.9° and an Io 
System III longitude angle of 48.6°, which are similar to the emission 4 observation conditions in Table II. Sodium was ejected uniformly from an 
assumed exobase of 2600 km radius with a velocity dispersion for a Maxwell— Boltzmann flux distribution, where in (a) — 1.3 km sec and 

</*, - x 10 8 atom cm 2 sec and in (b) v m = 2.0 km sec 1 and <A) = 1.8 X 10* atom cm' 2 sec" 1 (see text). 


two-dimensional sky-plane grid by extracting an average 
radial profile. This average radial profile (called the calcu- 
lated eclipse profile) will be denoted by the filled circles 
in Figs. 4-7. A one-dimensional east-west D 2 brightness 
profile (and also a corresponding column density profile) 
for the emission data and the sodium cloud image data is 
obtained from the larger two-dimensional sky-plane grid 
by selecting only the east-west grid elements that occur 
in the grid row containing Io. In Figs. 4-7, the calculated 
cast-west brightness and column density profiles are de- 
noted by filled triangles for the forward cloud profile and 
by filled squares for the trailing cloud profile. To construct 
an eclipse or east-west profile, monoenergetic model cal- 
culations are performed for 18 different nonuniformly 
spaced speeds ranging from 0.4 to 10 km sec" 1 . Profiles for 
speeds beyond 10 km sec" 1 are determined by an inverse 
speed extrapolation of the model results. The individual 
profiles for the different speeds are appropriately weighted 
for a given source flux speed distribution and then added 
to obtain the final profile. Model calculations are per- 
formed for an Io geocentric phase angle of 92.9° and an 
Io System III longitude angle of 48.6°. These satellite condi- 
tions arc similar to those for the emission 4 and eclipse 2 


observations of Table II, which are the observations closest 
to the eastern elongation point. This choice is also appro- 
priate for all the eclipse data, which has no discernible 
dependence on these two Io related angles, and for the Io 
sodium cloud image data which have east-west profile 
areas in Fig. 2 that are representative of the satellite near 
its orbital elongation points. Modeling analysis results are 
summarized in Table V and discussed below. 

For the first Maxwell-Boltzmann flux distribution in Fig. 
3 with a most probable speed of = 1.3 km sec" 1 (i.e., 
an exobase temperature of —1560 K) and with a flux <f > 0 
of 3.0 X 10 8 atoms cm' 2 sec" 1 (i.e., a total source of 
— 1.2 x 10 26 atoms sec” 1 ), the model calculated eclipse 
profile (filled circles) in Fig. 4a provides an excellent fit 
within the Lagrange sphere to the eclipse observations 
(open circles) and also compares very favorably with the 
east-west column density profiles calculated for the for- 
ward (filled triangle) and trailing cloud (filled squares). 
This fit verifies and is similar to the earlier 1 500 K Maxwell- 
Boltzmann flux distribution fit of Schneider et al (1991) 
noted in Section 2. Beyond the Lagrange sphere in Fig. 
4a, however, all three of these calculated profiles fall below 
the eclipse observations, which is considered less accurate 
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100.0 10.0 1.0 1.0 10.0 100.0 


RADIAL DISTANCE FROM IO (R, 0 ) 

FIG. 5. Model calculations for the east-west D 2 brightness profiles using a Maxwell-Boltzmann flux speed distribution. The east-west D? 
brightness profiles near lo in both the trailing and forward cloud directions as determined by the emission 4 data of Schneider et al (1991) arc 
shown by the open circles. The east-west profile envelopes in both the trailing and forward cloud directions as determined from the sodium cloud 
image data are shown by the shaded areas (see Fig. 2 caption). The descriptions for the calculated profile symbols, the sodium cloud mode! and 
plasma torus, and the Maxwell-Boltzmann flux distribution in (a) and (b) are the same as in the caption of Fig. 4. 
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FIG. 6. Model calculations for the eclipse data using a modified sputtering flux speed distribution. The atomic sodium column density profile 
near Io determined from the 1985 eclipse data by Schneider et al. (1991) is shown by the open circles. The model calculated column density profiles 
are shown by solid dots for the (cylindricatly averaged) corona, by solid triangles for the forward cloud along the east-west direction, and by solid 
squares for the trailing cloud along the east-west direction. These column density profiles were calculated using the Io sodium cloud model of 
Smyth and Combi (1988b) for their case C description of the plasma torus and for an Io geocentric phase angle of 92.9° and an lo System III 
longitude angle of 48.6°, which are similar to the emission 4 observation conditions in Table II. Sodium was ejected uniformly from an assumed 
exobase of 2600 km radius with a velocity dispersion for a modified sputtering flux distribution, where in (a) a - 3, v m = 1.0 km sec -1 , and <A,= 
3.2 X 10 K atom cm' 2 sec in (b) at - 7/3, = 0.5 km see \ and <A, = 4.2 x 10* atom cm" 2 sec -1 , and in (c) a = 2, v m = 0.4 km sec" 1 , and = 

4.7 x 10 H atom cm -2 sec -1 . 


1 ! 






IO S SODIUM CORONA AND SPATIALLY EXTENDED CLOUD 


71 


at these distances. At and beyond about 8 /? Io , the calcu- 
lated east-west forward (filled triangle) and trailing (filled 
squares) profiles rise above the calculated eclipse profile 
(filled circles) because the column density is no longer 
spherically symmetric about Io, with the forward cloud 
profile having the largest column density and showing a 
distinct change in its slope compared to the trailing cloud 
profile. The corresponding model profiles for the D 2 emis- 
sion brightness are given in Fig. 5a. For both the forward 
and trailing profiles, the calculated eclipse and calculated 
east-west profiles are in good agreement with each other 
inside the Lagrange radius, with a maximum brightness of 
about 200 kR near the exobase. The calculated east-west 
profile threads the three emission 4 data poinls for the 
forward cloud, but falls well below the emission 4 data 
points in the trailing cloud. For both the forward and trail- 
ing clouds at larger radial distances, the calculated east- 
west profiles fall well below the areas for both the forward 
and trailing cloud images. This behavior indicates that 
there is a large deficiency in the high-speed population for 
this source flux speed distribution. 

Model calculations were therefore performed for the 


second Maxwell-Boltzmann flux distribution in Fig. 3 with 
a higher most probable speed of u m = 2.0 km sec -1 (i.e., 
an exobase temperature of —3690 K) and with a flux 
<f ) o of 1.8 X 10 8 atoms cm 2 sec 1 (i.e., a total source of 
—0.75 X 10 26 atoms sec ') and are shown in Figure 4b and 
Figure 5b. For the D 2 emission brightness profiles in Figure 
5b, the calculated east-west profile now threads the center 
of the forward cloud image area for a radial distance up 
to about 70 R ]0 and the lower trailing cloud image area 
for a radial distance of about 25 R\ 0 before it fails off 
too steeply. This improved fit at larger radial distances, 
however, reduces the D 2 emission brightness at the exobase 
to about 80 kR in Figure 5b and causes the calculated 
eclipse profile in Figure 4b to fall below the measured 
eclipse profile for radial distances inside about 3 /? Io . The 
Maxwell-Boltzmann flux distribution therefore cannot fit 
both the corona profile near Io and the sodium cloud east- 
west profiles at large distances from the satellite. A flux 
distribution that has a broader dispersion with enhanced 
populations for both the low-speed and high-speed atoms 
is required. The three modified-sputtering flux distribu- 
tions in Fig. 3, which have a broader dispersion, are thus 




FIG. 7. Model calculations for the east-west D : brightness profiles using a modified sputtering flux speed distribution. The east-west D 2 
brightness profiles near Io in both the trailing and forward cloud directions determined from the emission 4 data of Schneider et ai (1991) are 
shown by the open circles. The cast- west profile envelopes determined from the sodium cloud image data are shown by the shaded areas (see 
caption of Fig. 2). The descriptions for the calculated profile symbols, the sodium cloud model and plasma torus, and the modified sputtering flux 
distribution in (a), (b). and (c) are the same as in the caption of Fig. 6. 
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TABLE V 

Summary Comparison of Modeled and Observed East-West Sodium Profiles for Different Flux Speed 

Distributions at Io’s Exobase 




Corona 

Forward Cloud 

Trailing Cloud (not null) 

Trailing Cloud (null) 

Observations: E/W Radial Interval (Ri 0 ): 

14-6 

6-10 

10-100 

6-10 

10-100 

6-10 

10-100 

Distribution 
Speed Peak 
(km/s) 

Exobase 
Source Rate 
(1026 aloms/s) 

Fits 

Eclipse 

Column 

Profile 

Fits 

Near Io 

Emission 

Profile 

Fits 
Far 
Cloud 
E/W Profile 

Fits 

Near lo 

Emission 

Profile 

Fits 
Far 
Cloud 
E/W Profile 

Fits 

Near lo 

Emission 

Profile 

Fits 
Far 
Cloud 
E/W Profile 

L MaxwellBoltzmannI 

r lux Distribution 








1.3 

1.24 

YES 

YES 

too low 

little low 

too low 

too low 

too low 

2.0 

0.75 

too low 

YES 

slightly low 

YES 

little low 

little low 

too low 


2, Colltsional Cascade Flux Distribution 
(classical sputtering) 

a = 3 10 1.32 YES YES tiny low YES tiny low little low too low 


(incomplete cascade: higher velocity tail) 


a - 7/3 0.5 

1.74 

YES 

YES 

YES 

YES 

YES 

little low 

too low 

a = 2 0.4 

1.90 

tiny low 

too high 

too high 

too high 

too high 

YES 

YES I 


considered in the remainder of the paper with model calcu- 
lations presented in Figs. 6 and 7. 

Model calculations for a classical sputtering flux distribu- 
tion (a = 3) and a modified-sputtering flux distribution 
(a = 7/3) are presented in Figs. 6a and 6b for the eclipse 
observations and in Figs. 7a and 7b for the east-west D 2 
emission brightness profiles. For these two flux distribu- 
tions, the most probable speeds are, respectively, 1.0 km 
sec -1 and 0.5 km sec" 1 , and the sodium fluxes <f > 0 are, respec- 
tively, 3.2 x 10 8 atoms cm" 2 sec -1 (i.e., a total source of 
— 1.3 x 10 26 atoms sec" 1 ) and 4.2 X 10 8 atoms cm' 2 sec" 1 
(i.e., a total source of —1.7 X 10 26 atoms sec" 1 ). From the 
exobase to radial distances of —8 R lo , just beyond the 
Lagrange radius, both sputtering flux distributions provide 
a very good fit in Figs. 6a and 6b to the observed eclipse 
column density profile (open circles) and correspond to 
an exobase D 2 emission brightness of about 150 kR in Figs. 
7a and 7b. For the classical sputtering flux distribution in 
Fig. 7a, the calculated D 2 emission brightness profile for 
the forward profile is slightly above the measured data 
point (open circles) inside the Lagrange radius, matches 
the two measured data points beyond the Lagrange radius, 
and then threads the forward cloud image area nicely be- 
tween about 20 R io and 80 R\ u before it falls too rapidly 
and drops below this area. An excellent fit for the forward 
profile is, however, provided by the modified sputtering 
distribution (a = 7/3) in Fig. 7b where the calculated D 2 
emission brightness profile matches the measured data 


points (open circles) both inside and beyond the Lagrange 
radius as well as nicely threading the forward cloud image 
area all the way to 100 R io . For the trailing cloud, the 
calculated D 2 emission brightness profile for the classical 
sputtering flux distribution in Fig. 7a matches the measured 
data point inside the Lagrange radius, is slightly below the 
two measured data points outside the Lagrange radius, 
and then threads the lower of the two trailing cloud image 
areas nicely between about 15 R lo and 35 R io before it falls 
too rapidly and drops below this area. An excellent fit for 
the trailing profile is, however, provided by the modified 
sputtering distribution (a = 7/3) in Fig. 7b where the calcu- 
lated D 2 emission brightness matches the measured data 
point inside the Lagrange radius, is slightly below the two 
measured data points outside the Lagrange radius, and 
then threads the lower (non-null) trailing cloud image area 
nicely all the way to 100 R ]0 . It is particularly noteworthy 
that the isotropic ejection of sodium from the exobase 
with a modified sputtering flux distribution with a = 7/3 
provides a complete fit to the combined eclipse, emission, 
and forward/trailing sodium cloud image profile data for 
this non-null condition from 1.4 to 100 /? [o . 

In order to fit the trailing cloud (upper area) profile for 
the directional feature at the null condition, it is then clear 
that a flux distribution is required with an even more en- 
hanced higher-speed population (—20 km sec" 1 ) than the 
modified sputtering flux distribution with a = 7/3. Since 
the modified sputtering flux distribution for a = 7/3 corre- 
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sponds to the limit of a single collision cascade process 
described by a Thomas-Fermi cross section (see Smyth 
and Combi 1988b), reducing the value of a to a smaller 
value becomes somewhat physically questionable but will 
be used here for the purposes of simply illustrating the 
impact of a more enhanced higher-speed sodium popula- 
tion in the model calculation. As discussed earlier, this 
higher-speed sodium is thought to be nonisotropically 
ejected from Io’s exobase and attributed to some combina- 
tion of direct collisional and lower-velocity charge ex- 
change ejection. Choosing the modified sputtering flux dis- 
tribution with a = 2 in Fig. 3 which has a most probable 
speed of 0.4 km sec 1 and selecting an isotropic exobase 
source rate of 1.9 x 10 26 atoms sec -1 (i.e., a flux of 
4.7 X 10 5 * * 8 atoms cm 2 sec" 1 ), the model-data comparison 
is shown in Fig. 6c for the eclipse column density and 
in Fig. 7c for the east-west D 2 emission brightness. The 
sputtering flux distribution provides a reasonably good fit 
to the observed column density data points in Fig. 6c with 
only a small departure very near the exobase and produces 
a column density profile beyond 10 R\ 0 that is significantly 
enhanced compared to the a = 7/3 case in Fig. 6b. In Fig. 
7c, this enhancement in the forward cloud is obvious, where 
the calculated D 2 emission brightness profile is significantly 
above the measured data points both inside and outside 
the Lagrange radius and is above or in the very top of the 
forward cloud image area all the way to 100 R Ur The 
additional enhanced high-speed population of the a = 
2 modified sputtering flux distribution is too large and 
therefore not consistent with the observed forward profile. 
In contrast for the trailing cloud in Fig. 7c, the calculated 
D 2 emission brightness profile matches the measured data 
points inside and outside of the Lagrange radius very well 
and then threads the upper of the two trailing cloud image 
areas nicely all the way to —90 R hy This demonstrates that 
the trailing cloud can be fitted with an enhanced higher- 
speed population of sodium atoms in the flux distribution. 
It also immediately demonstrates that the flux distribution 
at the exobase must be nonisotropic with the enhanced 
high-speed population weighted toward vector directions 
that will preferentially populate the trailing cloud. As dis- 
cussed in Section 2, this nonisotropic requirement for a 
flux distribution for speeds of —20 km sec 1 is similar to 
the conclusions reached by earlier modeling analyses (Pil- 
cher etai 1984; Smyth and Combi 1991; Wilson and Schnei- 
der 1995). 

5. DISCUSSION AND CONCLUSIONS * 

The composite spatial information for sodium obtained 
by combining the eclipse observations (radial distances 

from Ioof 1.4 to — 10 /? Io ), the emission observations (east- 

west distances of ±4 to ±30-40 /? lo ), and the sodium cloud 
observations (east-west distances of ±10 to ± 100 R lo ) has 



FIG. 8. Two-dimensional nature of the sodium column density near 
Io. Contours for the two-dimensional column density near Io are shown 
in the sky-plane of the earth as determined from the sodium cloud model 
calculation for the modified sputtering flux speed distribution described 
in Fig. 6(b) for at = 7/3. The vertical and horizontal directions are the 
projected directions that are, respectively, perpendicular and parallel to 
the semimajor axis of the Io’s orbital ellipse on the sky plane. Io's location 
and size are shown to scale by the black circle. The sodium column 
density contours in units of 10 11 atoms cm 2 are, from inside to outside, 
5, 2, 1, 0.3, 0.2, 0.1, and 0.05. 


been analyzed to extract a basic description for the flux 
speed distribution at the satellite’s exobase. An isotropic 
modified-sputtering flux speed distribution in Fig. 3 with 
a = 7/3, a most probable speed of 0.5 km sec" 1 , and a 
source strength of 1.7 x 10 26 atoms sec" 1 provides a very 
good fit to these composite observations when the direc- 
tional feature is either north or south and hence not con- 
tributing to the east-west profile of the trailing cloud. It 
is remarkable that these observations, acquired by a num- 
ber of ground-based programs over very different spatial 
scales and at different times during the 1976-1985 decade, 
are so self-consistent. Near Io, the two-dimensional sodium 
column density produced by this modified sputtering distri- 
bution as calculated by the sodium cloud model in the 
profile analysis above is shown in Fig. 8 and can be seen 
at larger distances from lo to become nonspherical and 
more confined near the satellite plane. This flattening near 
the satellite plane is the merging of the near Io corona 
into the sodium cloud and is caused naturally by orbital 
dynamics beyond the satellite Lagrange sphere where the 
gravity of Jupiter is dominant. The forward cloud portion 
of the east-west emission data profiles has a rather tightly 
confined slope that, in the absence of the trailing cloud 
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FIG. 9. Total source rate speed distribution function for sodium at Io’s exobase. The total source rate speed distribution function at Io’s exobase, 
in units of 10 26 atoms sec 1 (km/sec)" 1 , is composed of three separate source rate speed distributions as discussed in the text and is shown for two 
different source strengths for the higher-speed zenocorona source centered about 57 km sec -1 . The lower (solid line) and upper (dashed-dot line) 
curves correspond, respectively, to the sodium zenocorona higher-speed source rates of 2.2 X 10 26 atoms sec 1 (Smyth and Combi 1991) and 4 X 
10 26 atoms sec *. The two source rates for the higher-speed zenocorona source are shown to exhibit its typical time-variable source strength range 
of —2-4 x 10 2fi atoms sec -1 as reported by Flynn el ai (1994). The decomposition of the solid curve into its three separate source rate speed 
distributions is shown in the cutout and is determined by combining (1) the isotropic modified sputtering source rate distribution (dotted line in 
the cutout) for a = 7/3, v m = 0.5 km sec 1 and a source strength of 1.7 x 10 26 atom sec \ (2) the nonisotropic lower-speed source rate distribution 
(short dashed line in the cutout) for the sodium zenocorona and directional feature centered about 20 km sec \ with a source strength of 1.1 x 
10 26 atoms sec -1 as determined by Smyth and Combi (1991), and (3) the nonisotropic higher-speed source rate distribution (longer dashed line in 
the cutout) for the sodium zenocorona centered about 57 km sec -1 , with a charge exchange source strength of 2.2 x 10 26 atoms sec 1 as determined 
by Smyth and Combi (1991). 


enhancement at the null condition, is less steep and is 
brighter than the trailing cloud profiles. In order, however, 
to reproduce the extended east-west profile in the trailing 
sodium cloud when the directional feature is in the satellite 
plane (i.e., the null location), additional nonisotropic high- 
speed sodium is required. 

The sodium atoms ejected from Io’s exobase as de- 
scribed above by the modified sputtering flux distribution 
have speeds primarily in the range from 0 to a few tens 
of km sec 1 . This neutral flux distribution represents the 
spatially integrated effect of the incomplete collisional cas- 
cade process that occurs from the collisional interactions 


of heavy ions in the corotating plasma torus with neutrals 
in lo’s atmosphere. This flux speed distribution can be 
alternatively described as a source rate speed distribution 
by multiplying it by the satellite surface area. In addition to 
these ion-neutral elastic collisional encounters, resonance 
charge exchange between plasma torus sodium ions and 
neutral sodium in Io’s atmosphere is also responsible for 
producing a sodium source with higher speeds relative to 
Io. These speeds are centered about the corotational ion 
speed (—57 km sec' 1 ) relative to Io’s motion and have a 
dispersion reaching from several tens of km sec -1 to '-TOO 
km sec" 1 . Such high speed sodium (<80 km sec ! ) has 
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recently been observed near Io by Cremonese etui. (1992). 
As discussed in Section 2, this higher-speed nonisotropic 
source of sodium together with the lower speed (—15-20 
km sec 1 ) nonisotropic source for the directional feature 
form the source for the sodium zenocorona or magneto- 
nebula. Earlier modeling studies (Smyth and Combi 1991; 
Flynn et ai 1992) indicated that the higher-speed source 
was —2 X l() 26 atoms sec 1 while the lower speed source 
was ~1 X 10 2f> atoms sec -1 . More recent observations and 
analysis (Flynn et ai 1994) have shown that the higher- 
speed sodium source is time variable and in the range 
—2-4 X 10 26 atoms sec *. A total source rate speed distri- 
bution for sodium at Io's exobasc has hence been con- 
structed by combining the modilied sputtering source rate 
distribution determined in this paper with the two source 
rate distribution for the zenocorona as given by Smyth and 
Combi (1991) and is shown in Fig. 9. The lower (solid line) 
and upper (dashed-dot line) curves correspond, respec- 
tively, to the sodium zenocorona higher-speed source rates 
of 2.2 X H) 26 and 4 X 10 26 atoms sec" 1 . Total source rate 
speed distribution functions at Io’s exobase expected for 
other atomic species, such as K, O, and S, can be con- 
structed in a similar fashion by adopting the estimated 
source rates given by Smyth and Combi (1991). 

Future studies for the sodium flux speed distribution at 
Io’s exobase are anticipated using much larger data sets 
now available for east-west and north-south sodium emis- 
sion observations. It will then be possible to analyze the 
combined spatial and spectral information and refine the 
nonisotropic nature of the flux distribution and also to 
search for other possible east-west and System III related 
variations. Once this information is determined for sodium, 
the implications for the more abundant species in Io’s 
atmosphere will be particularly important in other related 
studies for the many-faceted and complex phenomena in 
the Io-Jupiter system. 
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Abstract. The first results for applying a three-dimensional multiscale ideal MHD model for the 
mass-loaded flow of Jupiter’s corotating magnetospheric plasma past Io are presented. The model 
is able to consider simultaneously physically realistic conditions for ion mass loading, ion-neutral 
drag, and intrinsic magnetic field in a full global calculation without imposing artificial 
dissipation. Io is modeled with an extended neutral atmosphere which loads the corotating 
plasma torus flow with mass, momentum, and energy. The governing equations are solved using 
adaptive mesh refinement on an unstructured Cartesian grid using an upwind scheme for MHD. 

For the work described in this paper we explored a range of models without an intrinsic magnetic 
field for Io. We compare our results with particle and field measurements made during the 
December 7, 1995, flyby of Io, as published by the Galileo Orbiter experiment teams. For two 
extreme cases of lower boundary conditions at Io, our model can quantitatively explain the 
variation of density along the spacecraft trajectory and can reproduce the general appearance of the 
variations of magnetic field and ion pressure and temperature. The net fresh ion mass-loading rates 
are in the range of -300-650 kg s -1 , and equivalent charge exchange mass-loading rates are in the 
range -540-1150 kg s' 1 in the vicinity of Io. 


1. Introduction 

Io’s interaction with Jupiter’s corotating inner 
magnetosphere has been studied for over 25 years. Io has a 
neutral atmosphere which is probably locally thick but rather 
uneven across its surface [Lellouch et ai, 1992; Baliester et 
ai, 1994]. The ultimate source for atmospheric gases appears 
to be the numerous active volcanoes on the surface. It is the 
energetic particle environment near Io which is responsible for 
the balance of the plasma heating, Joule heating, ionization, 
and surface and atmospheric sputtering, and which in some 
form drives the escape of the neutral atmosphere [Schneider et 
ai, 1989]. A realistic model for the flow of magnetospheric 
plasma, interacting with lo’s neutral atmosphere is necessary 
for a complete understanding of the global picture. 

The results of the measurements by the particle and field 
instruments on the Galileo Orbiter during the December 1995 
flyby of lo provide new and important information with which 
a realistic MHD simulation for the plasma interaction can be 
tested. Kivelson et al. [ 1 996a,b] have recently suggested that 
an intrinsic magnetic field for Io with a strength and direction 
which would nearly cancel the local Jovian Field near the 
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surface of Io, can explain many aspects of the Galileo 
magnetometer measurements taken during the December 1995 
flyby of lo. However, the large plasma densities and low 
temperatures measured at the same time by the Galileo plasma 
science package [Frank et ai, 1996] and the plasma wave 
instrument [Gurnett et ai, 1996], indicate that ion mass 
loading is significant and may provide an explanation for the 
magnetic field measurements without requiring an intrinsic 
field. 

Our understanding of the interaction of Io with Jupiter’s 
magnetosphere goes back to the discovery of lo-related 
decametric radio emission discovered by Bigg [1964] and the 
unipole inductor model which explains it [Goldreich and 
Lynden-Bell 1969]. There have been a number of theoretical 
studies of this interaction during the immediate post-Voyager 
era [e.g., Southwood et ai, 1980; Neubauer, 1980] (see also 
Hill et ai [1983] for an excellent post-Voyager review). 
Neubauer [1980] presented an analytical model of Alfven 
standing wave current system which connects current through 
the ionosphere of Io. Southwood et ai [1980] examined data 
from several Voyager instruments and examined the possible 
role of an intrinsic magnetic field for Io as a way to retain a 
robust enough ionosphere to provide enough conductivity for 
completing the Io-Jupiter current circuit. Early theoretical 
work was often done either in the context of a “thin” 
atmosphere [e.g., see Cloutier et ai, 1978] indicative of the 
surface temperature (-130 K), or “thick” extended neutral 
atmosphere [e.g., see Goertz, 1980] more indicative of volcanic 
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temperatures (-1000 K). Subsequent evidence [Sieveka and 
Johnson, 1985; Schneider et ai, 1991; Lellouch et ai . , 1992; 
Ballester et ai, 1994; Belton et ai, 1996] seems to indicate a 
mixed picture of the global atmosphere, which has a large 
extended corona like a thick atmosphere, but appears to be 

dominated by local major injection of hot (high speed) 

gas/dust plumes to high altitudes but only near active 

volcanic vents. Therefore, although the atmosphere is 
probably only locally thick, it still has a large extended 
neutral corona which might provide a sufficient source of 
impact ionization and photoionization. 

There have been a number of three-dimensional (3-D) 

numerical studies of the plasma flow past Io [Linker et ai, 
1988, 1989, 1991; Wolf-Gladrow et ai, 1987]. Because of the 
complexities in the problem and numerical difficulties, studies 
have been limited to restricted portions of the problem, have 
introduced fixes such as artificial viscosity, or have been run 
using less than fully physically realistic conditions. Linker et 
ai [1991] have performed the best example to date of 3-D 
MHD modeling oflo. In this model they used the two-step 
Lax-Wendroff method. These calculations were performed on a 
nonuniform but fixed structured grid, and artificial dissipative 
terms (artificial viscosity) had to be introduced to the basic 
MHD equations in order to make them computational stable. 

A new numerical method has been developed at the 
University of Michigan for solving the full 3-D MHD 
equations which resolves a number of difficulties that have 
plagued past attempts. The 3-D version of this new method, 
called multiscale adaptive upwind scheme for MHD (MAUS- 
MHD), has already been applied to the full comet-solar wind 
interaction problem and has been discussed in detail by 
Gombosi et ai [1996]. This interaction is quite similar to that 
expected for the Io-Jupiter-magnetosphere interaction in the 
absence of an intrinsic magnetic field for Io. A full Earth 
magnetosphere version of the same calculation has also 
recently been developed [De Zeeuw et ai, 1996], so the method 
is also able to treat intrinsic magnetic fields as well, if 
necessary. We can and have run test cases with an intrinsic 
field for Io, but in this paper, concentrate on trying to explain 
the existing data with nonintrinsic field models. 


2. Model Description 

In the 3-D MAUS-MHD model the calculation is done on 
an unstructured Cartesian grid using octree adaptive mesh 
refinement (AMR). This makes optimal use of computer 

resources allowing problems with disparate scales and 

unusual and varying geometries to be modeled globally, while 
still retaining sufficient resolution locally to capture shocks 
and other important features resulting from naturally occurring 
sharp gradients in the flow. The octree structure is a 
hierarchical cell structure, based on multiple generations of 
cubic parent cells which can be split into eight daughter cells. 
The adaptive data structure has been explained by DeZeeuw 
and Powell [1992]. The approximate eight- wave Riemann 
solver for ideal MHD is combined with a second-order 
MUSCL-type scheme [van Leer , 1979]. It uses a novel 
approach to handle the “div B problem” derived and 
discussed at length by Powell [1994] and Powell et ai [1995, 
also Multiscale adaptive upwind scheme for magneto- 


hydrodynamics; MAUS-MHD, manuscript in preparation, 
1998]. In this scheme div B/r is treated as a passive scalar so 
that any numerically generated nonzero contributions are 
convected away. Once div B = 0 is imposed as an initial 
condition to the problem it is satisfied to within truncation 
error. 

Our model of Io’s interaction with Jupiter’s magnetosphere 
is essentially similar to the comet problem of Gombosi et ai 
[1996]. The dimensionless conservative form of the ideal MHD 
equations can be written as 



where W and S are eight-dimensional state and source vectors, 
while F is an 8 x 3 dimensional flux diad. All quantities 
denoted by a tilde are normalized by physical quantities in the 
undisturbed upstream region 


t — aoot/Rj 0 

(2a) 

f= r/R Io 

(2b) 

p=p/p~ 

(2c) 

u = u/a„ 

(2d) 

P=p/P'*,&<~ 2 

(2c) 

^ v ^oP°° ac>0 ^ 

(2f) 


where t is time, r is radius vector, p is mass density, u is bulk 
flow velocity, p is pressure, Po ls the permeability of free space, 
and B is magnetic field vector. The subscript «> refers to values 
in the undisturbed upstream flow. Finally, Ri 0 is the radius of 
Io and aoo is the hydrodynamic sound speed in the undisturbed 
upstream plasma flow. 

The normalized state and flux vectors are written as 
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where e is the normalized internal energy density 
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The source vector in (1) corresponding to ion mass-loading, 
M, is of the same form as used by Gombosi et ai [1996], and is 
given in normalized units as 


ii 
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Qi is the ionization rate per unit volume evaluated at lo’s 
radius; Fj n is the ion-neutral collision momentum -energy 
transfer rate per unit volume, also evaluated at lo’s radius. In 
the simplest interpretation Qi(R lo /r) a is given by the neutral 
density, rin, divided by Tj, the ionization timescale for neutrals, 
and Fi n (Rio/r) a is given by n n kj n , where kj n is the ion-neutral 
momentum transfer collision rate, and n n is the neutral density. 
The choice for the exponent, a, is discussed below. For the Io 
cases presented here the net neutral flow velocity, u n , is taken 
to be zero. This could be generalized in the future to account 
for a net escaping neutral atmosphere at high altitudes. 

For the Io version of the MAUS-MHD model simulations 
were performed both for a pure conducting sphere with and 
without ion mass loading and for a sphere with fixed boundary 
conditions, as discussed in detail below. The simulation 
volume is 900 x 600 x 600 R lo , and there are 9 levels of 
refinement and 92,000 cells with sizes ranging from <0.1 Rj 0 
near Io to 50 Rio far upstream and downstream. We have 
chosen the upstream plasma conditions as follows from our 
best estimates of the Galileo December 1995 epoch, including 
the various recently published Galileo measurements (see 
Table 1). Since the calculation is performed in normalized 
units, there is some flexibility to rescale the results over a 
range of realistic conditions, but this cascades to other 
parameters [see, e g., Gombosi et ai, 1996]. We have 

investigated numerous versions of cases for a simple 
conducting sphere, a mass- loading ionosphere, and with an 
intrinsic magnetic field for Io. 

The formulation of mass loading for comets [Gombosi et ai, 
1996] permits us to incorporate both the effects of mass loading 


Table 1. Upstream Plasma Conditions 


Parameter 

Value 

Upstream plasma density 

3500 cm' J 

Upstream plasma temperature 

92 eV 

Upstream mean molecular mass 

22 amu 

Upstream magnetic field 

1800 nT 

Corotation flow speed 

56.8 km s' 1 

Alfvdnic Mach number 

0.4 

Mach number 

2.2 

Ratio of specific heats (y) 

1 667 


from new ionization (impact and photoionization) as well as 
momentum and energy loading from charge exchange and 
nonreactive ion-neutral collisions. In the comet calculation 
Gombosi et al. used the term “ion-neutral friction” to describe 
the effect of collisions which cause momentum and energy 
exchange but which result in no new ions. The dominant 
source of such an exchange of momentum and energy results 
from symmetric charge exchange reactions involving the 
dominant torus ions and atmospheric neutrals (0 + + O, S + + S, 
and S ++ + S). Linker et al [1991] appropriately included an 
ion mass loading rate and a “charge exchange rate” in their 
MHD models. However, many nonreactive collisions between 
ions and atmospheric neutrals also can exchange momentum 
and energy without yielding net new ion mass. These so- 
called knock-on collisions are important for producing the 
extended neutral corona of Io [Sieveka and Johnson , 1985; 

Johnson ; 1990]. Furthermore, it has been shown by Smyth and 
Combi [1997] that the distribution of sodium in the corona 
[, Schneider et ai, 1991] as well as the more familiar B-cloud is 
a natural outcome of the collisional cascade process 
(incomplete atmospheric sputtering) that is initiated by such 
collisions. For simplicity in the remainder of this paper this 
momentum-energy loading term is referred to as the charge 
exchange rate. 

On the basis of Voyager-era conditions, Smyth [1992] has 
estimated a total collision rate of -2 x 10 28 s‘ l between 
incident plasma torus ions and neutrals in lo’s atmosphere 
with about 1/3 yielding fresh ions deposited into the plasma 
flow and the rest contributing to momentum/energy exchange 
alone, i.e., the charge exchange rate. These estimates are based 
on extrapolations of knowledge gained about the extended 
sodium atmosphere of Io to the more abundant neutral species 
S, O, SO, and S0 2 . On the basis of their own estimates, Linker 
et al [1989, 1991] have explored MHD models with a 
comparable charge-exchange rate but with a fresh ion mass- 
loading rate nearly an order of magnitude smaller (~10 27 s' 1 ) 
than Smyth’s. Based on the then non-detection of optical/UV 
emissions close to Io, Shemansky [1980] provided an upper 
limit to the fresh ion mass-loading rate of ~10 27 s' 1 . Smyth and 
Shemansky [1983] later raised this to 4 x 10 27 s' 1 after 
subsequent detections of neutral O and S emissions from 
ground-based and space-based remote measurements. It is also 
likely that the torus is more robust in recent years, and that 
mass, momentum, and energy input rates may ail be higher now 
than 17 years ago. This appears to be borne out by 
comparisons between Galileo measurements during the 
December 1995 flyby and Voyager measurements [Frank et ai, 
1996; Gurnett et al ., 1996]. 

As is evidenced by the range of values which appear in the 
literature, making an independent estimate for ion mass loading 
and charge exchange rates is difficult, because they are based 
on a number of uncertain and indirect measurements and 
assumptions. With the availability of the Galileo PLS 
measurements [Frank et ai, 1996] of plasma density, 
temperature, and pressure we have adopted the different 
approach of using the net ion mass-loading and charge 
exchange rates as adjustable parameters and constraining the 
values based on comparison of our model results with the 
measurements along the spacecraft track. 

S0 2 appears to be the major primary constituent in lo’s 
atmosphere [Ballester et ai, 1994; Lellouch et ai, 1992], but 
observations still to date only provide hemispherically 
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averaged total column abundances. The sodium eclipse 

measurements of Schneider et ai [1991] are the only 
unambiguous information we have about the vertical 
distribution of neutral gas species in the corona of Io from 
about 1.4 to about 6 R lo . The interpretation of recent Hubble 
observations of various UV emissions of atomic S and O 
[Ballester et ai , 1996] is complicated by the fact that the 
spatial distributions of the emissions are determined by a 
combination of the mass distribution in the atmosphere with 
excitation by the plasma density and temperature environment. 
Theoretical and empirical models have explored the structure 
of the neutral atmosphere and its behavior to various gas 
sources and energetic heating and ejection processes, such as 
sublimation, volcanic ejection, surface and atmospheric 
sputtering, UV, plasma, and Joule heating, and LTE and non- 
LTE IR radiative cooling [Sieveka and Johnson , 1985; 

McGrath and Johnson, 1987; Moreno et ai, 1991; Strobel et 
ai, 1994; Pospiezalska and Johns ori, 1996; Wong and 
Johnson , 1995; Smyth and Combi, 1997]. Despite recent 
advances in observational and theoretical work, still only a 
rough sketch of the structure of the neutral atmosphere can be 
made. 

For our MHD simulations we have adopted ion mass, 
momentum, and energy loading source terms in the vicinity of 
lo as if proportional to the neutral density distribution which 
is scaled as a power law in distance (r/R Io ) with an exponent, a. 
We do not need to specify the neutral density itself, however, 
just the ionization mass and momentum/energy exchange rates 
and their dependence on distance. We have adopted the power 
law with an exponent of a - -3.5 as found by Schneider et ai 
[1991] for the trace species, sodium, outside about 1.4 R Io . 
Smyth and Combi [1997] have shown that this vertical 
distribution of sodium from the Schneider et ai [1991] eclipse 
observations and the more typical emission profiles of more 
spatially extended sodium can together be explained by a 
single modified incomplete cascade sputtering distribution, 
which results naturally from multiple neutral-neutral and ion- 
neutral collisions between torus ions and atmospheric 
neutrals. Therefore it is reasonable to expect a similar type 
distribution for neutral species in general. This power law is 
contrasted with the less steep inverse square-type power law 
distribution assumed in past MHD mass-loading calculations 
[Linker et ai, 1988, 1991]. For the mass-loading models 
presented here, we assume the neutral distribution has this 
dependence from an inner boundary at 150 km above the 
surface of Io, out to the Lagrange radius of 6 R io . 

If details about the ion composition and electron 
temperature in the torus, and the global spatial distribution 
and composition of the neutral atmosphere were available, it 
would be possible to make a better estimate of the spatial 
distribution of the ipn mass-loading rate and the charge- 
exchange rate (the collisional energy -momentum loading rate) 
needed to specify the source terms in the MHD equations. 
Typical cross sections [McGrath and Johnson, 1989] between 
oxygen and sulfur ions dominating the torus and O, S, and SO 2 
neutrals dominating the atmosphere show that typical 

symmetric charge exchange reactions have cross sections in the 
range of 20 to 40 A 2 whereas most impact ionization and 
elastic collisions [Johnson, 1990] have cross sections in the 
range of 1 to 20 A 2 . Electron impact of course only nets ion 
mass but is highly dependent on the electron temperature. 

The power law distribution may not be correct well inside 
of 1.4 R| 0 because the range of the observations of Schneider et 


ai [1991] started at this distance. At much lower altitudes the 
resulting ion mass-loading rate is a trade off between the likely 
rapid increasing of neutral density (probably more like an 
exponential scale height distribution), with the decreasing 
plasma flux (electron and ions) caused by recombination and 
plasma-neutral collisions themselves. Flaser et ai. [1997] 
have shown preliminary results of Galileo radio occultation 
measurements which indicate electron densities in excess of 
50,000 cm* 3 at the peak of the ionosphere, -50-80 km above the 
surface. This is similar to but higher than given by Kliore et 
ai [1974, 1975] from the Pioneer era, and again generally 
consistent with the more robust torus of the recent times. At 
this time, a fully coupled MHD-ionosphere calculation would 
be beset with many assumptions, for example, the neutral 
atmosphere density profile which is still uncertain to within a 
factor of a few. In any case the scopfc of this type of study is 
currently beyond the state of the art even for the Earth for 
which there is much more knowledge. 

For the models presented in this papier we have adopted two 
types of lower boundary conditions which represent opposite 
pictures of the type of conditions present above lo’s 
ionosphere. One is a perfectly conducting impenetrable 
sphere. In this case both the magnetic field and the plasma flow 
are reflected from the lower boundary. This represents an 
extreme case of not allowing either plasma flow or magnetic flux 
to penetrate the sphere, and we hereafter refer to this as a 
reflecting boundary condition. The second case corresponds 
to fixed boundary conditions where the magnetic field and 
plasma density and velocity are specified in ghost cells just 
below the boundary. In this case both magnetic flux and 
plasma flow can pass through the boundary. For the fixed 
boundary condition model presented here the plasma density 
is set to a value of 3 times the upstream torus density (or 
10,500 cm* 3 ) which is representative of the global spherical 
average indicated by radio occultation measurements at an 
altitude of 150 km. The plasma velocity is set to zero, being 
consistent with strong coupling to the neutral atmosphere. 
The choice of a value for the fixed magnetic field at the lower 
boundary requires some discussion. We hereafter refer to this 
case as fixed boundary conditions. 

Linker et ai [1991] presented model results which solved 
for the full MHD equations (mass, momentum, energy, and B 
field) outside of lo, and only the magnetic field with a constant 
conductivity inside the inner boundary, throughout the 
interior of Io. The internal conductivity was set in such a way 
to obtain a high enough integrated conductivity (10-150 
mhos) necessary for maintaining the Jupiter-Io current system. 
They point out that it is not really known where the current 
loops close, i.e., in the ionosphere but outside of Io, or even 
through the dense lower ionosphere and into a possible 
conducting interior. Linker et al. found, however, that the 
same exterior solutions could be obtained whether they 
performed the computationally intensive calculation of the 
magnetic field inside Io or whether they simply held the B field 
at the inner boundary constant at the background level. In fact, 
their full calculation showed that in any case “the magnetic 
field changed very little inside the sphere.” 

The physical interpretation for a constant magnetic field is 
that if Io does not have an intrinsic B field (i.e., from an internal 
dynamo), it would be expected that a value similar to the 
average imposed Jovian external field would be induced over 
time. We have estimated that the time-averaged value of the 
field at Io happens to be slightly less (-98%) than the 
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instantaneous value at Io’s position at the point of Galileo’s 
close approach. The same ratio was found either roughly using 
the variation from the old Goertz et al. [1976] magnetic field 
model or more carefully using the newest model of Connerney 
et al. [New model of Jupiter’s magnetic field constrained by 
the Io flux tube footprint, manuscript in preparation, 1998]. It 
is noteworthy that those models predict values that are overall 
too low and too high, respectively, compared with the actual 
values measured by Galileo [Kivelson et al., 1996b]. 

As mentioned already for this study we have left the ion 
mass-loading rate and the charge exchange rate (pure 
momentum/energy exchange) as adjustable parameters. We 
then undertook an iterative procedure to vary these quantities 
and compare our model results with the measured plasma 
density, pressure, and magnetic field along the Galileo 
spacecraft trajectory until a feasonable (albeit not perfect) 
match was obtained. Upstream we assume the aforementioned 
torus conditions. Downstream and to the sides of the 
simulation box we assume absorbing boundary conditions, i.e., 
mass, momentum, and energy, exit, and nothing comes back. 
Because of the adaptive grid technology we are able to place 
the upstream, downstream, and lateral boundaries far enough 
from Io so they do not affect on the solution. 

3* Model Results 

To serve as a baseline, we show results for a perfectly 
conducting sphere without mass loading, given realistic 
upstream pl££ma torus conditions. We refer to this case as no 
mass loading and reflective boundary conditions. This is 
shown mainly for comparison with the similar nonintrinsic 
field model of Kivelson et al. [1996b], and to serve as a point of 
comparison with the mass-loaded cases. For this case, the only 
parameters to be set are the upstream Mach number of 2.2 and 
the upstream Alv6nic Mach number of 0.4. In (5) the ion-mass 
and ion-neutral drag rates and F in are zero. The results are 
scaleable to other values of upstream magnetic field and 
conducting radius given constraints imposed by the Mach 
number (which does not have a great effect on the overall 
global structure of the flow) and the Alfvdnic Mach number 
(which is more critical). In any case our estimates of both the 
upstream Mach number and Alfvdnic Mach number are 
consistent with the measured Galileo flyby conditions and 
rigid corotation. 

We show results for two mass-loading cases where the 
particular values of the ion mass- loading rate and the charge- 
exchange rate serve as the adjustable parameters. We have run 
some models with larger exponents (-4 to -5), however, -3.5 is 
based on an actual measurements and is quite satisfactory. The 
goal is to be able to reasonably reproduce the variation of the 
plasma density, pressure (and temperature), and magnetic field 
along the spacecraft trajectory according to the results of the 
Galileo Io flyby results [Frank et al, 1996; Kivelson et al. , 
1996b; Gurnett et al, 1996]. The sensitive features in the 
plasma data to be reproduced are (1) the height of the plasma 
density peak along the center of the wake, (2) the reduction of 
the magnetic field, (3) total width of the disturbance in all 
quantities, (4) the reversal and side peaks of the pressure and 
temperature, and (5) the vector velocity magnitude and 
direction. Frank et al, [1996] have reported plasma velocities 
(44 km s" 1 ) that are below corotational speeds (56.8 km s" 1 ) 
even 4 Rj 0 from the central axis of the wake. However, these 


have a large enough uncertainty (-5 and +15 km s' 1 ) to be 
consistent with corotational speed, which seems more 
reasonable. 

Plates la- If show a comparison of slices of the magnetic 
field and flow streamlines in the planes parallel and 
perpendicular to the upstream magnetic field for the non-mass- 
loaded and both mass-loaded cases. In the non-mass-loaded 
reflective case one can see the strong influence of the Alfven 
wings in the parallel plane (la) [see Neubauer, 1980; Linker et 
al, 1991] as the flow diverts around the conducting sphere 
and as the disturbance is directed off at the expected angle (tan 
0 = l/M A ) where 0 is the deflection angle and M A is the 
Alfvdnic Mach number. In MHD the Alfvdn wings are the 
analog of vortex sheets of regular fluid dynamics [Kogan, 
1961]. In the parallel plane for the mass-loading case with 
reflective boundary conditions (lb), the effect of the inner 
boundary conducting sphere, which dominates the non-mass- 
loaded case, is weakened but still present. Finally, in the 
fixed-boundary mass-loading case (1c) the effect of the Alfven 
wings has all but disappeared. Plates Id- If show the 
comparable plots in the plane perpendicular to the upstream 
magnetic field. The December 1995 Galileo flyby trajectory lay 
very close to this plane. The diverted flow streamlines are not 
very different between any of the three cases, in fact for the two 
mass-loaded cases the flow fields in this plane are virtually 
identical. The small variations in the magnetic field magnitude 
are apparent. Note that in mass-loading fixed case (Plates 1c 
and If) the diversion of flow is similar in both planes. 

Plates 2a, 2b & 2c show the striking differences for the 
plasma ion density among the various cases. The dark lines 
show the outlines of the octree cell structure of the adaptively 
refined mesh. In the non-mass-loading reflective case (2a) the 
flow diversion around the sphere is evident. The density 
concentrates somewhat along the Alfv6n wings but the wake is 
essentially evacuated of material. The mass-loaded reflective 
case (2b) shows a large density pileup in the wake right 
behind Io and a fair dispersion of material spread in the 
direction parallel to the upstream magnetic field. This is 
caused by the strong influence of the reflective boundary 
conditions. The mass-loaded fixed boundary case (2c) also 
produces a plasma wake, but it is more concentrated along the 
central line of the wake, does not have the large extension 
along the field-parallel direction (as in Plate 2b), and is more 
extended along the wake. 

Direct comparisons of the model with the Galileo flyby 
data, along the trajectory are given in Figures la- Id. Model 
results are superimposed on the measurements for both mass- 
loading models with reflective and fixed boundary conditions. 
The reduction of the magnetic field in the wake along the 
spacecraft trajectory is shown in Figure la. Both of these 
models provide a similar breath and depth of the magnetic field 
disturbance as the pure intrinsic dipole model of Kivelson et 
al, [1996b]; i.e., the measured field disturbance is wider and 
deeper. The model with reflective boundary conditions yields 
a somewhat deeper central disturbance than the fixed boundary 
model, but the overall widths are nearly identical. 

Figure lb shows the plasma density compared with the 
measurements of Frank et al [1996]. Both the peak of the 
density at the center of the wake and the width of the 
disturbance are well reproduced by both mass-loading models. 
In fact, the height of the plasma density was one key factor in 
adjusting the mass-loading rate in our process of reproducing 
the measurements. The plasma pressure from the model is 
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compared with the Frank et al. data in Figure 1c and the 
temperature in Figure Id. The character of the pressure/ 
temperature disturbance is reasonably well reproduced. There 
is a rise to either side of a narrow central minimum. The depth 
and width of the central minima in pressure and temperature, 
which correspond well with narrow density peak, are well 
reproduced by both models. However, the measurements show 
a large increased temperature and pressure in the flanks of the 
wake, which roughly corresponds to the wider measured 
magnetic field disturbance, and which is not quantitatively 
reproduced by either mass-loading model. The models show 
smaller peaks which are closer to the central line of the wake, 
and the model with reflecting boundary conditions produces 
higher side peaks in both temperature and pressure. 

The two key attributes on which we focused fitting were the 
height of the density peak (as mentioned above) and the central 
depth of the minima in temperature and pressure. Although 
there is some interplay in adjusting both parameters, the 
density peak was more sensitive to changes in the mass- 
loading rate, while the depths of the pressure and temperature 
minima were more sensitive to the charge exchange rate. For 



y ( lo radii } 


the model results shown in Plates 1 and 2 and Figure 1, the 
total ion mass-loading rate corresponds to 0.84 x 10 28 and 1.8 
x 10 28 ions s' 1 for the fixed and reflective boundary condition 
models, respectively, where the mean molecular mass is 22 amu 
and the radial distribution varied as r' 3 5 . The charge exchange 
rate (due again to momentum/energy loading for which charge 
exchange is the dominant part) is equivalent to 1.5 x 10 28 and 
3.2 x 10 28 ions s" 1 , respectively, for comparison with the 
charge exchange rates given by Linker et al. [1991]. For the 
best match to the density peak and pressure reversal, the ion- 
neutral collision frequency was about two times the ionization 
rate. The 2 to 1 ratio in relative rates are in rough agreement 
with recent estimates by Smyth [1992]. 

This net ion mass loading corresponds to a fresh ion mass- 
loading rate of 300-650 kg s" 1 . From the momentum/energy 
exchange rate we can generate ah equivalent total apparent 
mass-loading rate which has contributions from both the true 
ionization rate as well as the momentum and energy loading 
from the combination of (mostly) symmetric charge exchange 
and non-reactive ion-neutral collisions of about 840-1800 kg 
s' 1 . This total loading rate is a factor of 2-4 larger than Smyth’s 




Figure 1. Comparison of the mass-loading models with Galileo particle and field measurements. The magnetic 
field measurements in (a) from Kivelson et al. [1996] are shown with the model values. The (b) plasma density, 
(c) pressure, and (d) temperature from the model are shown with those from Frank et al. [1996] obtained from 
the plasma ion measurements. The solid lines give the results for the model with mass-loading and reflective 
boundary conditions. The dashed lines give the results for the model with mass-loading and fixed boundary 
conditions. 
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No Mass Loading - Reflective Boundary 



Mass Loading - Fixed Boundary 



Plate 1 . Magnetic field and streamlines for non-mass-loading and mass-loading flows past Io. (a)-(c) The 
plane aligned along the upstream magnetic field direction. The draped field lines (generally vertical lines) and 
velocity stream lines (generally horizontal lines with arrows) are shown, (d)-(f) The plane perpendicular to 
the upstream magnetic field direction. As indicated, Plates la and Id correspond to the case with no mass- 
loading and reflective boundary conditions. Plates lb and le correspond to mass loading and reflective 
boundary conditions, and Plates lc and If correspond to the mass loading and fixed boundary conditions. 
The spectral color table indicates the magnitude of magnetic field. The strong influence of the Alfv6n wings is 
apparent on the non-mass-loaded flow (Plates la and Id), and is weakened with the addition of mass loading. 
The increased magnetic field perturbation in the wake is apparent in the mass-loaded cases. The spatial scale 
covers from -10 to +12 Io radii horizontally and to +8 Io radii vertically. 
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No Mass Loading • Reflective Boundary 



Mass Loading - Reflective Boundary 



Plate 2. Plasma density for the non-mass- loading and mass- 
loading flows past Io. (a) The plasma density for the models 
with no mass-loading and reflective boundary conditions, in 

(b) for mass-loading and reflective boundary conditions, and in 

(c) for mass-loading and fixed boundary conditions. The non- 
mass-loading model (Plate 2a) produces a rarefied region in the 
wake, contrary to the Galileo particle measurements. Both 
mass-loading models (Plate 2b and 2c) produce a strong 
density peak in the wake behind Io. The symmetry line 
through Io covers distances from -4 to +8 Io radii. 


Voyager-era estimates but may be consistent with the apparent 
overall larger Io atmosphere and torus activity in recent times 
compared with the Voyager era. 

The general flow pattern in both mass-loading models is 
similar to the typically suggested qualitative pictures. For 
example, see Schneider et at. [1991]. Very close to Io and 
especially in the extended wake region the flow is very slow (a 
few km s' 1 ) where the effects from the locally produced ions and 
from charge exchange are most important. Because the Galileo 
trajectory is nearly in the perpendicular plane and because 
there is not much difference between the two models in this 
plane a comparison of the velocity field with the Galileo 
results provides little discrimination between the two 
alternative boundary conditions. Figure 2 shows the flow 
velocity vectors along the spacecraft trajectory for both mass- 
loading models, which give virtually identical results. The 
flow directions and change in the velocity magnitude in the 
orbit plane are quite similar to those shown in Figure 1 of 
Frank et al. [1996]. Note again, however, that Frank et al. 
reported asymptotic velocity magnitudes far from the wake (4 
Io radii) somewhat below strict corotation (45 km s' 1 as 
opposed to 56 km s' 1 ). They warned that their uncertainties 
were asymmetric and large enough to include strict corotation, 
being dependent on the details of their assumed plasma ion 
composition. Our model results yield corotational velocities 
far from Io, with somewhat supercorotational velocities ( 1 3% 
enhancements) in the immediate flanks of the wake near Io, 
about 2 Io radii from the center of the wake. 

In ideal MHD the current density can be calculated as the 
curl of the magnetic field vector. Extracting quantities 
determined from components of field gradients near Io, where 
they are largest, is uncertain because of the presence of cut- 
cells (partial cubes which contain the inner boundary at one 
face), because of the conversion from cell-averaged parameters 
to node-values, and because of boundary effects. For the model 
runs shown here the cell structure was not optimized to obtain 
highly accurate currents. We can say at least qualitatively that 
the largest values of the current density in all three models 
occur near the body. As one might expect, currents associated 
with the Alfv6n wings are the strongest in the non-mass- 
loading reflective boundary model and weakest for the mass- 
loading fixed boundary model. In the mass-loading fixed 
boundary model we can also estimate the integrated current 
from the current density near Io to be on the order of a few x 10 6 
amperes. This is comparable to independently estimated values 
[Khurana et al., 1997] for both the Galileo and Voyager 
epochs. 

Preliminary pictures of the global shape of the local 
ionosphere of Io have been presented based on radio 
occultation measurements by Flasar et al [1997] and plasma 
wave measurements by Gurnett et al [1997]. This overall 
shape is dominated by the dense, stagnant, and heavily mass- 
loaded region indicated by the very low velocities which 
excludes the outer plasma torus flow, as seen in Plate 2. It is 
apparent then that the overall structure of Io’s upper 
ionosphere may be largely controlled by the combination of 
mass loading and plasma torus flow past Io. Furthermore, the 
upstream-dawn and downstream-dusk ionosphere profile 
asymmetry seen in the Pioneer 10 radio occultation 
measurements [Kliore et al, 1974, 1975], that are qualitatively 
similar to the preliminary Galileo measurement, is then most 
likely an upstream/downstream phenomenon rather than a 
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Figure 2. Velocity flow vectors along the December 1995 
Galileo flyby trajectory. Shown along the projected trajectory 
in Io’s orbit plane are the velocity vectors (direction and 
magnitude) from our mass- loaded MHD model calculations. 
The results for the fixed boundary model are shown but the 
reflective boundary model is virtually identical. Both the 
variations in direction and variation in magnitude of the 
velocity are quite similar to those by Frank et al. [1996] 
obtained from the Galileo plasma ion measurements. The 
absolute velocities from the measurements are uncertain at this 
time. 


dawn-dusk one. This has been noted previously by Kivelson 
et al. [1979]. Preliminary Hubble Space Telescope measure- 
ments of atomic oxygen and sulfur emissions near Io [Trauger 
et al, 1997; Ballester et al., 1997] also seem to show the 
presence of a bright wake. 

Our case of reflective boundary conditions considers the 
ionosphere to be a perfectly conducting and reflecting sphere 
at 150 km altitude above the surface, somewhat above the 
electron density peak [Flasar et al. 1997; Kliore et al. 1974 & 
1975], This is certainly an overestimate of even a large but 
finite conductivity. Our fixed boundary condition allows 
plasma and magnetic flux to pass, but may underestimate the 
effect of the underlying high ion densities (50,000 - 90,000 
cm' 3 ) at the ionospheric peak near 80 km where the ions are 
likely more coupled with the neutral atmosphere, increasing 
the conductivity. In any case these two limiting cases likely 
bracket the true influence of a more realistic lower ionosphere/ 
atmosphere on the outside flow. It is interesting to note that 
the choice of the boundary conditions does not drastically 
effect the character of the direct comparison with the Galileo 
measurements themselves, although it does substantially 
change the character of the flow above and below Io in the 
plane parallel to the upstream magnetic field where there are as 
yet no Galileo data. Because of the larger deflection in the 
field-aligned direction for the reflective compared with the 
fixed boundary conditions, a higher mass-loading rate is 
required for the reflective case. 

Note that we have already mentioned that our model with 


no intrinsic field, but with mass-loading and charge exchange 
can produce the magnetic field signature about as well as the 
vacuum superposition model of Kivelson et al [1996b]. 
Assuming almost the same background field composed of a 
superposition of Jupiter’s field and the intrinsic dipole as that 
of Kivelson et al [1996b], Khurana et al [1997] have recently 
presented a calculation of the magnetic field perturbation that 
would result from the plasma pressure as measured by Galileo 
PLS [Frank et al, 1996]. We refer here to the broad peaks in 
plasma pressure (and temperature) which appear in the flanks at 
about 2 R Io from the center of the wake. There is no 
corresponding density increase at this location, but it does 
coincide with the broad magnetic field signature. They 
conclude that given the validity of their assumptions only 
-30% of the total field perturbation can be due to plasma effects 
as separate from an intrinsic magnetic field. In contrast both 
our models without an intrinsic magnetic field are able to 
account for roughly 2/3 of the total disturbance in magnetic 
and all of the density peak. 

Although our mass-loading profile does account for the 
density peak, and our charge exchange profile accounts for the 
temperature and pressure minimum in the center of the wake 
they do not produce the large broad temperature and pressure 
peaks and the corresponding larger magnetic field disturbance 
in the flanks of the wake. If as Khurana et al. [1997] found, the 
enhanced plasma pressure in the flanks of the wake can account 
for the general breadth and depth of the magnetic field 
disturbance (remember this is just outside the density peak), 
then it is equally reasonable to conclude that if we could add a 
process or mechanism which would produce the enhanced 
temperature (and resulting pressure in absence of increased 
density), it should similarly also produce the correct enhanced 
magnetic field disturbance in our model without an intrinsic 
magnetic field, just as it does for their model. 

Clearly then what is required is to alter our formulation to 
include a process which adds energy (or possibly subtracts 
less energy) in the flanks of the wake and/or upstream of the 
flanks. We consider two immediately apparent possible 
candidates: (1) the spatial dependence of the mass-loading and 
charge exchange rates, and (2) the formulation of the ion pickup 
process. In absence of including a detailed ion-chemical 
network, we adopted mass-loading and charge exchange rates 
which are both eveiywhere proportional to the measured 
neutral sodium radial profile of Schneider et al. [1991]. If the 
mass-loading rate and charge exchange rates are in reality not 
proportional to one another everywhere, it might be possible 
to enhance the temperature in the flanks of the wake. This 
increased temperature would necessarily have a corresponding 
increased pressure (where there is no increased density) and by 
the same argument as Khurana et al [1997] there should also 
be a corresponding further decreased magnetic field. It is worth 
noting that our models produce somewhat super-corotational 
velocities enhanced by 13% just at 2 Ri 0 from the center line of 
the wake, exactly where the temperature is enhanced. A 
different velocity dependence on collision cross sections 
would in general be expected for the different relevant 
ionization and charge exchange (and knock-on) processes. 

The precise forms for the ion-mass-loading and ion-neutral 
momentum -energy exchange (charge exchange) as functions of 
location around Io are in reality set by the interplay of a whole 
series of ion, neutral, and electron impact reactions which 
depend on the detailed composition of both the neutral 
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atmosphere and plasma. Including all the effects of this kind of 
detailed picture for the Io flow problems has already been done 
for the comet solar-wind flow problem in terms of successful 
comparisons with in situ spacecraft data [Gombosi et ai, 
1996] and 2-D ground-based spectra and images for H 2 0 + 
[Haberli et ai, 1997]. In the comet problem, although much 
more complicated in terms of total species and reactions than 
for Io, the details of the neutral atmosphere structure and 
composition are also much better understood than for Io. 

The second possibility would be to reexamine the premise 
behind the formulation of the pick-up process. As mentioned 
in the model description, the ion pick-up process adopted 
assumes ions are produced initially in a pick-up ring and 
shortly thereafter scattered into a shell distribution by wave- 
particle interactions. Recent Hubble Space Telescope images 
of the Jupiter aurora with the new STIS instrument [Clarke et 
ai, 1997] show that the footprint of Io is accompanied by a 
trailing footprint of the wake itself extending up to one-quarter 
to one-half way around Jupiter. This implies that ions picked 
up at Io, are still being pitch angle scattered into the loss cone 
for 3 to 6 hours, and still precipitating along wake-crossing 
flux tubes into Jupiter’s ionosphere for a quarter to a half a 
Jupiter rotation. Those picked up very close to Io, where 
densities are high, would be scattered immediately. This could 
also explain why the models provide an excellent match to all 
the measurements in the densest regions. The inclusion of 
pickup ions remaining in ring distributions would 
substantially alter the formulation of the MHD equations. A 
multispecies model that could account for a mixed ring and 
shell populations is beyond our current capability. It is 
certainly well beyond any modeling that has yet been 
published. It is possible that the inclusion of this process (W. 
Smyth, private communication, 1997) might have desired effect 
on the magnetic field disturbance. 


4. Summary 

We have demonstrated the application of an ideal MAUS- 
MHD model to the interaction of Jupiter’s corotating plasma 
torus with Io’s atmosphere. The model solves the single-fluid 
MHD flow in 3-D on an adaptively refined mesh. Because of 
its numerical robustness and the ability to adapt the grid 
structure to the flow, the heavily mass-momentum-energy 
loaded flow in the supersonic and sub-Alfv6nic regime 
relevant for the interaction of Jupiter’s plasma torus has been 
solved without the need to introduce artificial dissipative 
terms. 

We find that we can reproduce many of the observed 
structures of the plasma wake behind Io as measured by various 
Galileo instruments. We obtain a good quantitative match to 
the plasma density peak, the plasma pressure and temperature 
minima at the center of the wake, and the vector velocity plasma 
flow directions .and magnitude variations. We obtain a 
qualitative match to the overall pressure and temperature 
profiles. We reproduce the magnetic field disturbance about as 
well as the intrinsic dipole magnetic field model of Kivelson et 
ai [1996b]. The distant flanks of the disturbance indicate a 
more decreased magnetic field as well as corresponding 
increased plasma pressure and temperature which appear to be 
coincident and possible correlated. Although the responsible 
mechanism is not yet known, it may be due to decreased 


cooling in the flanks from velocity-dependent collision cross 
sections (which are not yet taken into account) or from the 
influence of pickup ions that are still in ring distributions. 
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